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I.     INTRODUCTION 

The  dynamic  response  of  a  submarine  under  casualty  conditions  constitutes  a 
crucial,  and  frequently  limiting,  factor  in  establishing  the  vehicle's  submerged  operating 
envelope.  As  basic  casualty  conditions  in  the  context  of  this  work,  we  mean  the  loss  of 
control  surface  and/or  propulsion  system  response,  and  a  flooding  casualty  where  the 
boat  must  either  be  brought  to  the  surface  or  stabilized  to  a  new  operating  depth.  Of 
particular  significance  is  the  study  of  the  ability  of  the  boat  to  recover  from  a  control 
surface  jam.  There  exist  several  factors  which  determine  the  severity  of  such  a  situation 
as  well  as  the  recovery  procedures;  the  initial  conditions  (speed,  depth,  pitch  angle,  etc.) 
during  the  jam,  the  actual  control  surface  angles,  reversing  time  and  backing  power  of 
the  propulsion  system,  the  ability  to  blow  ballast,  and  the  time  between  recognition  of 
casualty  and  initiation  of  proper  recovery  procedures.  To  this  end,  it  is  crucial  that  we 
have  a  clear  understanding  of  the  dynamics  of  the  boat  during  an  emergency  situation. 

The  traditional  methods  for  establishing  dynamic  stability  of  motion  concentrate 
mainly  on  eigenvalue  analysis  during  small  perturbations  around  nominal  straight  line 
paths  Clayton  and  Bishop  [Ref.  1].  Two  indices  are  utilized,  a  stability  index  G^  for  the 
vertical  plane  and  G^  for  the  horizontal  plane  Roddy  [Ref.  2].  In  terms  of  the  slow 
motion  derivatives,  these  indices  are  given  by: 
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Positive  values  for  these  indices,  whose  usage  dates  back  to  the  1950's,  indicates  motion 
stability  in  the  corresponding  plane.  The  underlying  assumption  in  these  criteria  is  that 
during  normal  straight  line  motions,  the  coupling  between  horizontal  and  vertical  plane 
motions  is  relatively  weak,  and  can  be  neglected.  Although  vortex  shedding  and  flow 
separation  introduces  a  certain  degree  of  coupling,  the  above  assumption  has  been  proven 
quite  useful  in  design  and  analysis.  However,  for  a  high-speed  fast-maneuvering 
submarine  operating  at  the  extremes  of  her  submerged  operating  envelope,  or  during 
emergency  situations,  the  above  assumption  of  uncoupled  motions  breaks  down.  High 
amplitude  motions  may  take  place  in  all  six  degrees  of  freedom,  and  the  nonlinear 
interactions  between  the  various  modes  of  motion  become  more  pronounced.  Therefore, 
we  have  to  carefully  consider  the  motion  characteristics  allowing  for  coupling  between 
horizontal  and  vertical  planes. 

The  implications  of  nonlinear  effects  and  coupling  are  numerous.  In  the  case  of  roll 
motion,  which  is  one  of  the  most  critical  responses,  there  is  growing  evidence  of 
complicated  dynamics  and  chaotic  response  under  certain  excitations  Falzarano,  Shaw  and 
Troesch  [Ref.  3],  Taz  Ul  Mulk  and  Falzarano  [Ref.  4]. 


The  motion  of  a  submarine  in  the  ocean  environment  involves  some  of  the  most 
complex  fluid-structure  interaction  problems.  For  example,  the  steady  and  unsteady 
characteristics  of  the  stem,  and  especially  with  regard  to  the  effect  of  the  wake  and 
rolled-up  body  vortices  on  propulsor  unsteadiness,  and  consequences  for  noise  and 
hull  vibrations  are  important  problems  of  great  concern.  At  present,  too  little  is 
known  in  depth  about  the  flow  at  the  stem  of  ships  and  submarines  and  the 
footprints  of  their  ensuing  wakes  in  homogenous  and  stratified  media  Sarpkaya 
[Ref.  5]. 

Typical  vortex  stmcture  about  a  submerged  body  can  be  seen  in  Figure  1 ,  which  is  taken 

from  Lugt  [Ref.  6]. 


TSTTTT 


Figure  1.     Typical  Vortex  Stmcture  About  a  Submerged  Body  (From  Lugt,  1981). 

When  an  axisymmetric  body  moves  at  a  sufficiently  large  angle  of  attack,  the 
boundary  layer  vorticity  may  lift  off  the  suri'ace  and  sheets  of  vorticity  are 
convected  downstream  while  rolling  up  into  streamwise  vortices  on  the  leeward 


side  of  the  body.  This  separation  and  roll-up  phenomenon,  known  as  the  crossplane 
separation,  occurs  on  almost  all  maneuvering  bodies  at  sufficiently  high  angles  of 
attack.  The  roll-up  of  the  vortices  and  their  subsequent  evolution  are  extremely 
important  for  the  proper  design  of  a  submerged  body  and  more  important  for  the 
determination  of  the  interaction  of  the  vortices  with  the  stem  and  the  propulsion 
system.  It  is  this  interaction  that  determines  the  maneuverability  and  control  of  a 
submerged  body  Sarpkaya  [Ref.  5]. 

In  the  case  of  submarine  motions,  there  is  evidence  of  bifurcation  phenomena  and 
extreme  sensitivity  of  response  to  initial  conditions  and  control  actions  during  emergency 
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Figure  2.    Typical  Simulation  Results. 


ascent  scenarios  such  as  recovery  from  a  dive  plane  jam.  As  motivation  for  the  analysis 
that  follows  we  present  typical  simulation  results  in  Figure  2.  The  time  simulation  is  in 
terms  of  vehicle  roll  angle  $  versus  time  for  2  %  excess  buoyancy  with  the  buoyancy 
force  located  1  %  of  the  vehicle  length  forward  of  the  center  of  gravity,  -6  degrees  dive 
plane  angle,  and  0. 1  feet  metacentric  height.  Four  different  recovery  actions  are  shown, 
all  parametrized  by  the  applied  rudder  angle  in  degrees.  It  can  be  seen  that  although 
zero  rudder  angle  appears  to  bring  the  vehicle  roll  angle  back  to  zero  in  the  shortest 
time,  the  response  does  not  persist.  Small  nonzero  values  for  the  rudder  angle  develop 
excessive  roll  angles  in  a  divergent  way,  while  larger  rudder  angles  reduce  the  amount 
of  roll.  Clearly  the  vehicle  effective  rudder  angle  is  largely  affected  by  the  amount  of 
vorticity  and  currents  in  the  flow  field,  and  its  actual  value  can  not  be  known  exactly. 
Since  situations  like  the  one  presented  in  the  graph  should  be  avoided,  we  need  to 
develop  a  mechanism  for  assessing  those  regions  in  the  parameter  space  where  the 
response  can  not  be  simulated  with  confidence. 

Stability  in  emergency  ascent  has  been  studied  in  the  vertical  plane.  In  Booth  [Ref. 
7]  and  [Ref.  8],  the  vehicle  response  was  distinguished  into  either  a  nearly  vertical  ascent 
or  predominantly  forward  motion.  In  Papoulias  and  McKinley  [Ref.  9],  it  was  found  that 
the  above  distinction  is  not  always  meaningful  as  a  result  of  the  many  parameters  that 
affect  the  problem.  This  latter  study  maintained  some  vertical  plane  restrictions,  although 
it  indicated  the  potential  existence  of  pitchfork  bifurcations  which  led  to  coupled  out-of- 
plane  solutions.  In  this  work,  we  relax  the  requirement  for  vertical  plane  motions,  and 
we  analyze  the  stability  properties  of  all  possible  steady  states  in  six  degrees  of  freedom. 


We  employ  a  combination  of  singularity  theory  Golubitsky  and  Schaeffer  [Ref.  10], 
bifurcation  theory  Guckenheimer  and  Hohnes  [Ref.  11],  and  numerical  continuation 
methods  Seydel  [Ref.  12]  in  order  to  capture  all  steady  state  solutions  that  are  physically 
admitted  by  the  coupled  nonlinear  equations  of  motion.  The  primary  bifurcation 
parameters  used  are:  amount  and  location  of  excess  buoyancy,  diveplane  and  rudder 
deflection,  and  the  metacentric  height.  Solution  branching  is  shown  to  occur  in  various 
forms,  including  single  and  multiply  connected  pitchfork  bifurcations,  separation  of 
solutions,  hysteresis,  and  teardrop  branches.  Dynamic  loss  of  stability  in  the  form  of 
Poincare-Andronov-Hopf  bifurcations  to  periodic  solutions  is  also  identified.  It  is  shown 
that  for  certain  ranges  of  parameters,  these  periodic  solutions  persist  in  the  vicinity  of 
inverted  pendulum  steady  states.  Finally,  we  summarize  our  results  in  the  form  of 
bifurcation  graphs  which  identify  parameter  regions  with  qualitatively  different 
asymptotic  response  characteristics.  For  demonstration  purposes,  all  computations  in  this 
work  are  performed  for  the  Swimmer  Delivery  Vehicle,  a  17.4  feet  vehicle,  for  which 
a  complete  set  of  hydrodynamic  and  geometric  properties  is  available  Smith,  Crane  and 
Summey  [Ref.  13].  Unless  otherwise  specified,  all  results  in  this  work  are  presented  in 
dimensional  form,  linear  dimensions  in  feet,  velocities  in  ft/sec,  angular  deflections  in 
degrees,  and  time  in  seconds. 


n.     PROBLEM  FORMULATION 

A.      PROBLEM  STATEMENT 

Controlling  emergency  ascent  situations  on  submersible  vehicles  such  as  dive  plane 
jam  recovery  is  of  concern  to  the  submariners.  In  order  to  control  such  situations,  one 
must  first  be  able  to  predict  the  dynamic  response  of  positively  buoyant  submersibles. 

Dynamic  response  equations  of  motion  describe  the  maneuvering  characteristics  of 
submersible  vehicles  for  six  degrees  of  freedom.  These  equations  assume  constant 
coefficients  for  hydrodynamic  forces  and  moments  approximated  by  zero  frequency  added 
mass  and  damping  terms  plus  the  quadratic  terms  for  drag  forces.  The  constant 
coefficients  vary  for  each  vehicle  and  depend  on  such  things  as  vehicle  body  shape, 
location  and  magnitude  of  vehicle  buoyancy,  position  of  bow  and  stem  planes,  position 
of  rudder,  vehicle  speed,  vehicle  mass  characteristics,  vehicle  hydrodynamic  coefficients, 
propeller  ipm  and  control  surface  inputs.  This  thesis  uses  the  equations  of  motion  and 
hydrodynamic  coefficients  for  a  submerged  Mark  IX  Swimmer  Delivery  vehicle  (SDV) 
developed  by  Smith,  Crane  and  Summey  [Ref.  13]  to  forecast  the  dynamic  behavior  of 
submersibles  in  a  positively  buoyant  condition. 

For  submarines  where  high  amplitude  motions  may  take  place  in  all  six  degrees  of 
freedom,  nonlinear  interactions  between  the  various  modes  of  motion  may  become  more 
pronounced.  In  particular,  there  is  growing  evidence  of  bifurcation  phenomena  during 
high  speed  maneuvering  and  emergency  ascent  scenarios  such  as  recovery  from  a  dive 


plane  jam.  In  these  cases,  it  is  no  longer  true  that  decoupled  linear  analysis  techniques 
are  sufficient  and  one  is  forced  to  consider  the  true  character  of  six  degrees  of  freedom 
motions.  So,  in  this  thesis,  six  degrees  of  freedom  equations  of  motions  have  been  used 
to  investigate  branching  behaviors  of  the  equations  of  motions.  This  is  done  under  the 
condition  of  changing  certain  system  parameters  while  keeping  others  constant. 

B.      EQUATIONS  OF  MOTION 

The  equations  of  motion  given  below  are  referenced  to  a  right-hand  orthogonal  axis 
system  fixed  in  the  body  as  shown  in  Figure  3.  The  origin  of  this  body  axis  system  is 
located  in  the  vehicle's  X-Z  plane  of  symmetry  on  the  mid-body  centerline  behind  the 
vehicle's  nose.  Positive  directions  for  control  surface  deflections,  forces  and  angular 
moments,  and  linear  and  angular  velocity  components  are  shown  in  Figure  3.  In  the 
following  equations,  differentiation  with  respect  to  time  is  denoted  by  a  dot  over  a 
quantity.  The  six  degrees  of  freedom  equations  of  motion  for  a  submarine  in  surge, 
sway, heave,  roll,  pitch,  and  yaw,  respectively,  are: 

inx[u-vT+wg-X(,(g^+T^)  -^-yaipq-r)  +Zg(pr+g)  ]  ^X^+X^f+X^       (1) 

mx  [v+ur-wp+X(;{pq+r)  -Vcip^-^-T^)  +Zg(gr-p)  ]  =Y^+7^+y(-.   (2) 

mx  [w-uq+vp+X(jipr-q)  +yG(<?r+p)  -z^^ip^+q^)  ]  =Z„+Zy,-*-Zc        (3) 
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+mx[xQ(y+ur-wp)  -yfJiu  -vr +h'^)]  ^Njj+N^'^N^ 
In  these  equations,  the  left  hand  sides  represent  inertial  forces  and  moments 
(Newton's  Law)  and  the  right  hand  sides  model  the  external  forces.  Subscript  H  reflects 
hydrodynamic  contributions,  W  buoyancy  and  weight  effects,  C  forces  arising  from 
control  surface  (rudders,  dive  planes,  and  bow  planes)  actions,  and  the  rest  of  the 
symbols  are  based  on  standard  notation  and  will  be  explained  at  the  end  of  this  section. 
The  hydrodynamic  radiation  and  viscous  forces  are  expressed  as: 

--P  r'^iCj,  h(x)(v^xrf+C„  bix)iw-xqf]^^^^^dx  ^'^ 
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These  are  given  in  customary  form  of  series  expansion  in  terms  of  the 
hydrodynamic  coefficients  and  cross  flow  integral  terms  which  are  integrated  over  the 
entire  length  of  the  body  and  represent  quadratic  forces.  The  cross  flow  velocity  U^f  is: 

U^=[(v+xr)^H'^-xqf]^^  (13) 

Hydrodynamic  restoring  forces  and  moments  are  due  to  the  vehicle  weight  W  and 
buoyancy  B,  and  are  given  by: 

Xyy=-(W-B)sinB  (14) 

Y^=(W-B)cosQsiik^  (15) 

Z^=iW-B)cosecos^  (16) 

A^,^=(ycH^-y^)cosecos<|)  -(z^  W'-z^)cosesin<|)  (17) 
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M^=  -(XfjW-XgB)cosQcos^  -(Zg  W^-z^)sme  (18) 

N^=(Xf.W-XgB)cosQsm<^  +Cy^Fr-)'^)sine  (19) 

Forces  and  moments,  due  to  control  surface  deflections,  are  reflected  as  added  drag 
in  surge,  while  in  sway,  heave,  pitch,  and  yaw  they  are  directly  proportional  to  the 
control  surface  deflections. 


(20) 


Y,-Y,u\  (21) 


Zc=u\Z,b^^Z,6^  (22) 


Kc=0  (23) 


Mc=u\M,b^^M,b^  (24) 

Nc=N^u^6^  (25) 

Usually,  control  surface  deflections  are  kept  intentionally  small,  and  the  linearity 
assumption  in  Equations  (21),  (22),  (24) ,  and  (25)  remains  valid.  Unlike  the  surface  ship 
case,  the  roll  moment  K^  is  zero  for  a  submersible  since  the  rudder  is  centered  with  the 
vehicle  center  plane.  The  hydrodynamic  coefficients  in  equations  (7)  through  (12)  are 
functions  of  the  frequency  of  motion,  or  what  amounts  to  the  same  thing  functions  of 
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the  maneuver  at  hand.  In  this  work  slowly  varying  reference  motions  are  studied  and 
therefore,  it  is  assumed  that  they  remain  constant  and  equal  to  their  zero  frequency  limit. 
It  should  be  emphasized  though,  that  the  constant  coefficient  assumption  would  break 
down  in  studies  related  to  the  fast  motions  under  the  action  of  first  order  wave  forces  in 
the  case  of  a  nearly  surfaced  submarine.  Another  important  assumption  in  this  study  is 
that  they  are  assumed  to  be  constant  throughout  the  range  of  vehicle  angle  of  attack. 
Ordinary  maneuvering  models  are  usually  validated  for  angles  of  attack  ±  15  degrees. 
For  higher  angles  of  attack,  the  cross  flow  drag  terms  Coy  and  Cj,z  dominate  the 
response,  and  they  are  functions  of  the  side  and  angle  of  attack.  Considering  them  to  be 
constant  does  not  alter  significantly  the  behavioral  characteristics  and  qualitative 
bifurcation  results  that  are  derived.  Similarly  ,  Cdy  and  Cd2  are  functions  of  speed  due 
to  the  Reynolds  number  effect  on  cross  flow  drag.  This  is  more  pronounced  in  small  size 
unmanned  untethered  vehicles,  whereas  for  submarines  the  cross  flow  drag  terms  remain 
relatively  constant  over  the  entire  speed  range. 

Since  the  equations  of  motions  are  referred  to  an  axis  system  that  is  fixed  in  the 
vehicle,  and  thus  translates  and  rotates  with  it,  the  orientation  and  the  position  of  the 
moving  body  axis  system  relative  to  fixed  inertial  reference  system  must  be  specified. 
The  orientation  of  the  body  axis  system  with  respect  to  the  inertial  reference  system  is 
defined  by  the  standard  Euler  Angles  ^(yaw),  0(pitch),  and  *(roll).  The  rotation 
sequence  from  the  inertial  reference  system  to  the  body  axis  system  is  ^,  9,  €>  as  shown 
in  Figure  4.  So  to  complete  the  model,  we  need  the  expressions  for  the  Euler  Angle 
rates  of  change; 
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(1)  Vehicle-Centered 
Gravity-Directed  System 
para  1  lei  to  inertia  1 
axis  system. 


(2)  Horizontal  Flight 
Reference  System  derived 
from  XqYqZo  by  rotation  a- 
bout  Z  through  yaw  angle  ^. 


y. (yMj 


(3)    Roll-Resolved   Flight 
Reference   System  derived 
from  X"Y"2"   by   rotation 
about  Y"   through  pitch 
angle  9. 


(k)    Vehicle  Body  Axis   Ref- 
erence  System  derived    from 
X'Y'Z'    by   rotation  about   X' 
through   roll   angle  ^. 


Figure  4.     Unit  Sphere  Development  of  Euler  Angles 
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4)=  p+  q  sin(|)  tane  +  r  cos(t)  tanO  '^"^ 

6=^  cos4)-r  sin4)  '^'^ 

^=  q^^  +  r^^  (28) 

cos6         cos6 

Finally,  it  is  assumed  that  propulsion  is  inoperative  and  the  propeller  is  rotating 
freely.  Therefore,  propulsive  forces  are  not  included  in  Equations  (1)  through  (6).  The 
driving  mechanism  for  the  vehicle  is  its  excess  buoyancy,  B-W  >  0.  The  problem  is  then 
to  assess  the  asymptotic  dynamic  characteristics  of  the  system  during  this  condition  of 
free  positive  buoyancy  ascent. 

The  model  presented  above  can  be  written  in  its  state  space  form  by  selecting  as 
state  variables; 

X^—U    y  ^— V    ,  x^—w    ,  x^—p  f2Q^ 

where  the  first  six  describe  the  system  motion  and  the  last  two  its  geometry.  Notice  that 
the  yaw  angle  "^  does  not  affect  the  equations  of  motion  and  is,  therefore,  not  included 
in  (29).  Angle  i'  can  be  computed  from  (28),  once  the  time  histories  of  the  state 
variables  have  been  obtained. 

In  a  compact  vector  form  the  state  equations  can  be  written  as, 

i    =  /(x)  (30) 

where  bold  face  indicates  a  vector. 
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Major  variables  and  parameters  as  defined  by  Smith,  Crane,  and  Summey 
[Ref.  13]  are  given  below: 

1.  Dynamic  Variables 

u,v,w  -  Linear  velocity  components  of  vehicle  with 

respect  to  orgin  of  body  axes  system  relative  to    fluid. 

p,q,r  -Angular  velocity  components  of  vehicle  with 

respect  to  body  axes  system  relative  to  inertial 
reference  system. 

X,Y,Z  -  Hydrodynamic  force  components  along  body 

axes. 

K,M,N  -Hydrodynamic  moment  components  along  body 

axes. 

2.  Mass  Distribution  Parameters 

m  -  Mass  of  the  flooded  vehicle,  including  the  mass  of  the 

entrained  fluid. 

W  -  Weight  of  the  flooded  vehicle  ,  including  the  weight  of 

the  entrained  fluid.  (=gm  ; where  g  is  the 
acceleration  of  gravity). 

V  -  Displacement  volume  of  the  vehicle. 

B  -  Buoyancy  force  acting  on  the  vehicle  (pgV).  This  is 

independent  of  the  inertial  mass  distribution  of  the 
submersible  vehicle,  including  whether  or  not  it  is 
flooded  . 

Xo  ,  Yq  ,  Zq      -   Coordinates  of  the  CG  (center  of  gravity)  in  the  body 
axis  system  (Figure  3).  These  will  depend  on  the  mass 
distribution  of  the  vehicle,  including  the  mass  of  the 
entrained  fluid. 

Xb,  Yb,  Zb         -  Coordinates  of  the  (center  of  buoyancy)  in  the  body  axis 
system  (Figure  3).  These  are  independent  of  the  mass 
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distribution  system,   but  may  vary  with  the  addition  or 
removal  of  external  appendages. 

I,,  ly,  I2  -  Moments  of  inertia  about  the  body  system  axes, 

including  the  entrained  fluid. 

Ixy>  Ik>  lyz  "  Products  of  incitia  about  the  body  system  axes, 

including  the  entrained  fluid. 

3.      Remaining  Parameters 

p  -  Mass  density  of  fluid  medium. 

b(x),  h(x)  -  Width  and  height  of  vehicle  in  its  xy  and  xz 

planes,  respectively,  at  location  x  measured  in  the 
body  axes  system  (Figure  3).  These  quantities  are 
required  in  the  integral  defming  the  cross  flow 
forces  and  moments  in  the  equations  of  motion. 

Xnose>  x^  "  Coordinates  of  vehicle  nose  and  tail  as 

measured  in  body  axis  system.  (Figure  3). 

8,  A>^r  ■  Stemplane,  bowplane  and  rudder  deflection 

angles  in  radians  (Figure  3). 

C.      STEADY  STATE  CONDITIONS 

Steady  state  conditions  are  achieved  when  the  submersible  reaches  constant  linear 
and  angular  velocities  which  must  assume  fmite  values.  Therefore,  the  body  fixed  linear 
accelerations  and  body  fixed  angular  accelerations  will  be  zero.  Likewise,  the  vehicle 
will  have  reached  constant  angles  of  roll  <f>  and  pitch  6  ,  making  the  derivatives  of  them 
equal  to  zero.  When  these  values  are  put  into  Equations  (1)  through  (6),  and  substituting 

(j)  =  e  =  0 

in  Equations  (26)  through  (28)  the  steady  state  values  of  the  angular  velocities  can  be 
found  as: 
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p  =-^  sin8  (31) 

q  =  ^  siiMJ)  cos6  (32) 

r  =  i|f  cos4>  cos  6  (33) 

Combining  Equations  (31),  (32),  (33)  with  the  equations  (1)  through  (6)  yields  a 
system  of  nine  coupled  nonlinear  algebraic  equations  in  the  form 

fix)  =0  (34) 

where  the  overbar  denotes  an  equilibrium  solution. 

Equation  system  (34)  can  be  solved  for  the  steady  state  values  of  u,  v,  w,  p,  g,  r, 
(/>,  6,\J/.  This  is  a  highly  nonlinear  system  of  equations,  and  it  may  exhibit  solution 
branching  and/or  multiple  solutions  McKinley  [Ref.  14].  Here  i^  is  not  zero  at  steady 
state  because  taking  ^=0  at  steady  state  will  restrict  the  analysis  to  the  vertical  plane. 
In  this  analysis,  motion  in  all  six  degrees  of  freedom  is  studied. 

Equation  system  (34)  can  be  written  in  the  following  form 

/(J,X)=0  (35) 

Here,  X  denotes  any  one  of  the  system  parameters.  By  system  parameter,  we  mean  the 
effects  of  control  surfaces,  propeller  rpm  the  values  of  buoyancy  and  weight,  etc. 
Therefore,  our  system  of  Equations  (35)  represent  a  multiparameter  problem.  This 
multiparameter  problem  can  be  solved  by  keeping  all  parameters  fixed  except  one 
(denoted  by  X).     This  system  of  equations  for  various  values  of  X  can  be  solved 
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continuation  method.  The  continuation  method  will  be  explained  in  the  following 
section. 

D.      PRINCIPLES  OF  CONTINUATION 

The  system  of  nonlinear  'algebraic'  equations 

f(y,X)=0  (36) 

serves  as  a  basis  for  the  discussion.  Here  y  denotes  an  n-dimensional  vector.  No 
generality  is  lost  by  restricting  attention  to  Equation  (36).  Both  steady  state  solutions  of 
ODE's  and  PDE's  are  solved  by  approximating  them  by  such  system  of  nonlinear 
equations. 

Assume  that  at  least  one  solution  of  Equation  (36)  has  been  calculated.  Let  us 
denote  this  first  solution  by  (y'  ,  X,).The  continuation  problem  is  to  calculate  further 
solutions, 

(y' ,  ^2) ,  (y^  X3) ' (^'^ 

until  one  reaches  a  target  point,  say  X  =  X^.  (The  superscripts  are  not  to  be  confused 
with  exponents.) 

The  j-th  continuation  step  starts  from  (an  approximation  of)  a  solution  (y",  Xj)  of 
equation  (36)  and  attempts  to  calculate  the  solution  (y*"^',  Xj+,)  for  the  next  X,  namely 
Xj+i- 

With  predictor-corrector  methods,  the  step  j-*j  + 1  is  split  into  two  steps: 
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Predictor  step 


Corrector  step 


see  Figure  3  Seydel  [Ref.  12]. 


[y]* 


^^^"•Vi' 


>(y^^'.^^^,) 


(y^Aj 


Figure  5.    Schematic  Showing  Predictor  Corrector  Method  (From  Seydel,  1988) 

In  this,  thesis  a  predictor  corrector  method  is  used.  In  general,  the  predictor  is  not 
a  solution  of  equation  (36).  The  predictor  merely  provides  an  initial  guess  for  corrector 
iterations  that  home  in  on  a  solution  of  equation  (36).  The  major  portion  of  the  work  is 
either  on  the  predictor  step  (resulting  in  an  approximation  close  to  the  branch)  or  on  the 
corrector  step  (if  the  predictor  has  produced  a  guess  far  from  the  branch).  The  distance 
between  the  two  consecutive  solutions  (y,  Xj)  and  (y*',  \j+,)  is  called  the  step  length  or 
step  size  Seydel  [Ref.  12]. 

It  is  obvious  that  first  of  all,  we  need  a  solution  procedure  for  the  particular  class 
of  problems  under  investigation.  In  this  work,  a  Newton-like  method  has  been  used  to 
fmd  the  solutions  of  the  problem.  Therefore,  calculation  of  the  first  solution  is  often  the 
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most  difficult  part.  Initial  guesses  may  be  based  on  solutions  to  simplified  versions  of 
the  underlying  equations.  Sometimes  choosing  special  parameters  (often  equal  to  zero) 
leads  to  equations  that  are  easily  solved.  In  complicated  cases,  a  hierarchy  of 
simplifications  must  be  gradually  relaxed  stqp  by  step,  the  solution  of  each  equation 
serving  as  an  initial  guess  to  the  following  more  complicated  equation.  Solving  a 
sequence  of  equations  with  diminishing  degree  of  simplification  until  finally  the  full 
equation  solved  is  called  homotopy  Seydel  [Ref.  12], 

Equation  system  (35)  can  be  simplified  by  assuming  the  angle  of  yaw  (^),  its 
derivative,  and  p,  q,  r  to  be  zero.  By  doing  this,  the  motion  of  the  vehicle  will  be 
restricted  to  the  vertical  plane.  In  this  case,  the  resulting  equations  can  be  solved 
analytically.  These  analytical  solutions  were  done  in  McKinley  [Ref.  14].  In  this  work 
these  solutions  and  certain  numerical  simulation  results  have  been  used  as  initial 
conditions  to  the  continuation  runs. 

The  continuation  approach  does  not  guarantee  that  all  possible  solutions  can  be 
located  in  a  specific  example.  In  particular,  changes  to  detect  isolated  branches  are 
small.  Therefore,  in  this  work  there  might  be  other  solutions  which  had  not  been  located 
by  the  continuation  algorithm  which  is  used.  In  this  thesis,  a  continuation  method  is 
implemented  BIFPACK  [Ref.  15].  A  careful  combination  of  analytical  techniques, 
continuation  methods,  and  numerical  simulations  are  used  to  ensure  that  all  steady  state 
solutions  are  captured. 
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m.     CONTINUATION  IN  6s 

A.      GENERAL 

In  this  chapter,  the  results  of  continuation  in  terms  of  the  stemplane  angle  will  be 
discussed.  As  was  mentioned  before,  it  is  very  important  to  have  an  initial  condition 
before  starting  a  continuation  study.  In  his  master's  thesis,  McKinley  solved  the  system 
of  equations  in  the  vertical  plane  [Ref.  14].  In  that  work,  he  found  that  when  the 
longitudinal  center  of  buoyancy  (Xgb)  was  equal  to  -1  percent  of  the  vehicle  length,  for 
a  certain  value  of  stemplane,  deflection  the  real  part  of  one  of  the  system  eigenvalues 
became  zero  while  the  others  stayed  negative.  By  using  this  fact,  it  was  assumed  that 
there  is  a  pitchfork  bifurcation  in  our  system  for  the  above  specific  values  of  system 
parameters.  In  this  work,  the  results  of  his  thesis  were  used  as  initial  conditions  in  the 
continuation  runs  to  find  the  vertical  plane  solutions.  Since  these  results  were  the 
analytical  solutions  of  the  simplified  system  equations  by  restricting  the  motion  in  vertical 
plane,  we  were  able  in  this  work  to  validate  the  continuation  runs.  In  continuation  runs, 
the  system  of  equations  Equation  (35)  was  solved  for  changing  values  of  the  primary 
parameter  X  (in  this  case  X=6s).  After  a  few  continuation  runs,  a  pitchfork  bifurcation 
was  found.   This  pitchfork  bifurcation  is  explained  in  the  following  section. 
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B.      PITCHFORK  BIFURCATION 

Suppose  we  have  the  nonlinear  system  of  state  equations, 

X  =  f{x)  (38) 

we  know  that  the  equilibrium  points,  x  of  the  system  are  defined  by: 

f(x)=0  (39) 

This  is  a  nonlinear  system  of  algebraic  equations,  and  it  may  have  multiple 
solutions  in  x,  which  means  that  the  nonlinear  system  may  have  more  than  one  position 
of  static  equilibrium.  If  we  pick  one  equilibrium  x  ,  we  can  establish  its  stability 
properties  by  linearization.  The  linearized  system  becomes 

x=Ax  , 

where  A  is  the  Jacobian  Matrix  of  f(x)  evaluated  at  x, 

^   ax'- 

and  the  state  x  has  been  redefined  to  designate  small  deviations  from  the  equilibrium  x, 

x^x-x. 

As  long  as  all  eigenvalues  of  A  have  negative  real  parts,  we  know  that  the  linear  system 
will  be  stable.  This  means  that  the  equilibrium  x  will  be  stable  for  the  nonlinear  system 
as  well. 

The  question  we  ask  ourselves  next  is,  what  happens  if  one  real  eigenvalue  of  the 
linearized  matrix  A  is  zero?  The  interesting  case  here  is  when  the  rest  of  the  eigenvalues 


23 


have  all  negative  real  parts,  otherwise  x  is  unstable  and  the  problem  is  solved.  If  the 
case  of  a  zero  eigenvalue  appears  to  be  too  specialized  to  be  of  any  practical  use, 
consider  this:  Assume  that  f(x)  depends  on  one  physical  parameter,  and  that  physical 
parameter  is  allowed  to  vary  over  some  range.  Then  it  is  clear  that  A  will  depend  on 
that  parameter  and  as  the  parameter  varies,  it  is  possible  that  one  real  eigenvalue  of  A 
will  become  zero  for  a  specific  value  of  the  parameter.  Our  problem  is  then  to  establish 
the  dynamics  of  the  nonlinear  system  as  one  real  eigenvalue  of  A  crosses  zero;  i.e. ,  goes 
from  negative  to  positive.  As  the  solutions  evolve  with  time,  things  are  interesting  only 
along  the  direction  of  the  eigenvector  that  corresponds  to  the  critical  eigenvalue  (the  one 
that  crosses  zero).  Along  the  rest  of  the  directions  in  the  state  space,  everything  should 
converge  back  to  equilibrium;  remember  that  we  assumed  that  all  remaining  eigenvalues 
of  A  have  negative  real  parts. 

We  can  see  then  that  it  is  possible  to  approximate  our  original  system  by  a  one 
dimensional  system,  which  is  much  easier  to  analyze.  The  dynamics  of  the  two  systems 
will  be  qualitatively  similar.  The  formalization  of  the  above  reduction  procedure 
constitutes  what  is  known  as  center  manifold  reduction,  or  normal  form  computation  in 
nonlinear  analysis.  So,  let  us  see  what  happens  for  the  case  of  a  zero  eigenvalue  by 
using  a  (typical)  first  order  system, 

x=Xx-x^     , 

for  X  >  0  there  are  two  nontrivial  equilibria,  x=  ±  VX.  The  transition  of  stability  is 
illustrated  by  Figure  6.   We  can  summarize  the  stability  as  follows: 
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•  For  X  <  0  only  the  trivial  equilibrium  exists  and  is  stable. 

•  For  X  >  0  the  trivial  equilibrium  becomes  unstable  and  a  pair  of  symmetric  stable 
equilibria  are  generated. 


stable 


stable 


'   unstable 


table 


Figure  6.    Supercritical  Pitchfork  Bifurcation 


This  phenomenon,  the  loss  of  stability  of  an  equilibrium  and  the  generation  of  additional 
equilibrium  states,  is  called  a  pitchfork  bifurcation  and  is  very  common  in  nature;  Euler 
buckling  of  a  beam  is  a  very  typical  example.  In  particular,  the  above  case  is  referred 
as  the  supercritical  pitchfork,  this  is  a  rather  benign  loss  of  stability  since  upon  loss  of 
stability  of  the  trivial  equilibrium  the  additional  nearby  equilibrium  states  are  stable. 
Occasionally,  the  above  case  is  referred  to  as  a  soft  loss  of  stability,  since  for  small 
values  of  X  beyond  its  critical  value,  the  final  steady  state  of  the  system  does  not  differ 
much  from  the  nominal  (trivial)  steady  state. 
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As  a  second  example,  consider  a  "similar"  system  as  before,  the  linear  part  remains 
the  same,  and  the  nonlinear  part  x^  changes  sign, 

equilibria  and  stability  behavior  for  this  equation  are  shown  in  Figure  7.  There  is  again 
a  loss  of  stability  at  the  bifurcation  point  (y,  X)  =  (0,  0),  but  in  contrast  to  the  first 
example  of  Figure  6  there  is  no  exchange  of  stability.  Instead,  the  stability  is  lost  locally 
at  the  bifurcation  point.   We  can  summarize  the  stability  as  follows  : 

•  For  X  >  0  only  the  trivial  equilibrium  exists  and  is  unstable. 

•  For    X   <  0  the  trivial  equilibrium  becomes  stable,  and  a  pair  of  symm.etric 
unstable  equilibria  are  generated. 


unstable 


table   ? 


N 


/ 


unstable 


unstable 


Figure  7.     Subcritical  Pitchfork  Bifurcation 

This  case,  which  is  also  shown  in  the  figure,  is  called  a  subcritical  pitchfork.  A 
comparison  with  the  first  case  reveals  that  this  is  a  much  more  serious  loss  of  stability 
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case.    Upon  loss  of  stability  of  the  trivial  equilibrium  position,  there  is  no  other  stable 
equilibrium. 

Now  let  us  perturb  the  second  example,  considering  the  equation 

Xx  +  x3-hP=0  (40) 

An  analysis  reveals  that  in  the  perturbed  case  the  bifurcation  is  destroyed. 
Differentiating  equation  (40)  with  respect  to  X  gives: 

Since  x=0  is  not  a  solution  of  equation  (40),  there  is  no  horizontal  tangent 
dx/dX=0  to  the  solution  curve  x(X).  Checking  for  a  vertical  tangent  means  considering 
the  dependence  of  X  on  x  and  diifeientiating  with  respect  to  x.  This  yields 

-^=0        foi       X=-3x2  , 
dx 

which,  substituted  equation  (40),  in  turn  yields  the  loci 

(xa)  =  [(|)^  .-3(|)^] 
of  the  vertical  tangents.   There  are  three  solutions  of  equation  (40)  for 

X     <     -3(|)^ 

and  the  branching  diagram  looks  like  Figure  8  (solid  lines).   The  point  with  a  vertical 
tangent  is  a  turning  point.  For  /S  <  0,  the  curves  in  Figure  8  are  reflected  in  the  X-axis. 
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As  jS  -*  0,  the  curves  approach  the  pitchfork  (dashed  curve).  For  an  arbitrary 
perturbation  /3,  the  occurrence  of  a  turning  point  is  typical.  Bifurcation  is  the  rare 
exception  because  it  occurs  only  for  i8=0  . 


(/S   >   0) 


/ 


Figure  8.    Perturbed  Pitchfork  Bifurcation 

In  our  case,  we  found  a  pitchfork  bifurcation.  We  found  that  the  trivial  solutions 
of  our  system  come  from  in-plane  (vertical  plane)  solutions.  After  the  pitchfork 
bifurcation,  we  had  a  branching  in  the  steady  state  solutions.  The  physical  meaning  of 
these  branches  is  that  they  correspond  to  out  of  plane  motion  of  the  vehicle  at  the  steady 
state.  To  get  these  branches,  we  used  the  numerical  simulation  results  as  an  initial  guess 
for  the  continuation  run.  Figure  9  shows  a  typical  pitchfork  bifurcation  that  we  found  is 
a  supercritical  pitchfork  bifurcation.  In  this  and  all  subsequent  runs  in  this  chapter,  5s  is 
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our  main  bifurcation  parameter  (X)  and  5r  is  the  perturbation  parameter  (/S).  In  Figure 
9  and  all  subsequent  figures,  solid  lines  correspond  to  the  symmetric  system  (6r=0 
degrees)  dashed  lines  to  6r=0.5  degrees,  dotted  lines  to  6r  =  2.5  degrees,  and  dash-dot 
lines  to  5r=  5  degrees.  As  we  can  see  in  this  figure,  it  illustrates  all  the  typical 
characteristics  of  a  pitchfork  bifurcation. 


TYPICAL  PITCHFORK  BirURCATION  FROM  CONTINUATION  in   DS 
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Figure  9.        Pitchfork  Bifurcation  From  Continuation  Runs 


C.      CONDITIONS 

1.       Defining  Additional  Terms 

a.    Excess  Buoyancy  BB 

Excess  buoyancy  is  defmed  as  5B  =  B  -  W  where  B  is  the  submersible' s 
total  buoyancy  and  W  is  the  submersible' s  total  weight 
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b.  Longitudinal  Center  of  Buoyancy  ,  Xcb 

The  longitudinal  center  of  buoyancy  is  defmed  as  Xgb  =  Xq  -  Xr  where 
Xg  is  the  longitudinal  center  of  gravity  with  respect  to  the  body  fixed  axis  and  Xg  is  the 
longitudinal  center  of  buoyancy  with  respect  to  the  body  fixed  axis. 

c.  Vertical  Center  of  Buoyancy,  Zcb 

The  vertical  center  of  buoyancy  is  defined  as  Zgb  =  ^g  -  ^Cg  where  z^  is 
the  vertical  center  of  gravity,  with  respect  to  the  body  fixed  axis,  and  Zg  is  the  vertical 
center  of  buoyancy  with  respect  to  the  body  fixed  axis. 

2.      Assumed  Conditions 

a.  Lateral  Centers  of  Gravity,  yc,  and  Buoyancy,  yg 

The  lateral  center  of  gravity  and  center  of  buoyancy  are  assumed  to  be 
on  the  same  plane(yG  =  ya  =  0). 

b.  Propeller  Speed  ,  n  (revolution  per  minute) 
The  propeller  speed  is  assumed  to  be  zero  (n=0). 

c.  Propeller  Coefficients  ,  Kj,^  and  N^^ 

From  Smith,  Crane,  and  Summey  [Ref.  13], the  propeller  coefficients  are 
zero(Kp,^  =  N^^  =  0). 

d.  Vertical  Center  of  Buoyancy,  Zcb 

The  vertical  center  of  buoyancy  is  assumed  to  be  positive. 
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D.      STEADY  STATE  RESULTS 

Figures  10  through  18  show  steady  state  solutions  for  the  eight  state  variables, 
namely;  u,  v,  w,  p,  q,  r,  i,6.  In  Figures  10  through  18,  solid  lines  correspond  to  6r  = 
0  degrees,  dashed  lines  correspond  to  dr  =0.5  degree,  dotted  lines  correspond  to  6r  = 
2.5  degrees,  and  dash-dot  lines  correspond  to  6r  =  5  degrees.  Figures  16  and  17  show 
steady  state  solutions  for  roll  angle  in  two  different  axis  scales.  In  all  these  cases, 
stemplane  action  is  continuation  parameter  (i.e.  ,X)  and  dr  is  the  perturbation  parameter 
(i.e.  ,|3).  Values  of  the  other  parameters  are  kept  fixed:excess  buoyancy  6B  =  2  %  of 
the  vehicle  weight  (W); deflection  of  bow  planes,  Sb  =  0; location  of  horizontal  and 
vertical  centers  of  buoyancy,  Xg  =  Zg  =  0;  location  of  vertical  center  of  gravity  ,Zgb  = 
0.1  feet;and  location  of  longitudinal  center  of  buoyancy;  Xob=  -1  %  o  the  vehicle  length. 
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Figure  10.     Steady  State  Surge  Velocity 
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SWAY  VELOCITY  V  vs  OS 
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Figure  11.     Steady  State  Sway  Velocity 
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Figure  12.     Steady  State  Heave  Velocity 
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Figure  13.  Steady  State  Roll  Velocity 
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Figure  14.    Steady  State  Pitch  Angular  Velocity 
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YAW  ANGULAR  VELOCITY  r  va  DS 
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Figure  15.     Steady  State  Yaw  Angular  Velocity 
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Figure  16.     Steady  State  Roll  Angle 


15 


20 


1    . 
1    . 
1    . 

j- \ 

:   /   ;  -N  :     i                    ^<-T-'  : 

••"''      ^-'''y           /    -i            • 

.1         .  t 

.1 

<         .1 

.1 
/ .         .1 

15 


20 


34 


ROLL  ANGLE  PHI  vs  DS 
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Figure  17.      Steady  State  Roll  Angle 
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Figure  18.     Steady  State  Pitch  Angle 
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It  can  be  seen  that  all  horizontal  plane  variables,  v,  p,  r  and  *,  exhibit  the  typical 
characteristics  of  pitchfork  bifurcation.  For  5r=  0,  there  exists  a  dive  plane  angle  6s 
between  -5  and  -10  degrees,  such  that  the  in  plane  solution  becomes  unstable  and 
symmetric  stable  out  of  plane  solutions  appear.  Nonzero  values  of  6r  result,  naturally, 
in  out  of  plane  solutions  for  the  entire  range  of  8s.  Multiple  solutions  appear  also  at  a 
certain  value  of  6s,  called  turning  point,  which  is  a  function  of  6r. 

The  vertical  plane  variables  u,  w,  q,  and  6  are  even  functions  of  horizontal  plane 
motions  and,  as  a  result,  they  exhibit  "solution  separation"  of  the  in  plane  motions  for 
values  of  6s  lower  than  the  bifurcation  point.  Of  special  interest  is  an  examination  of  the 
steady  state  values  of  the  pitch  angle  $  and  roll  angle  *,  since  these  are  directly  related 
to  the  vehicle  orientation.  A  clear  examination  of  the  possible  steady  state  solution,  so 
is  the  pair  (x-^,  x-*).  When  6r=  0,  it  can  be  seen  from  Figure  18  that  the  stable  in- 
plane  solution  exhibits  a  steady  state  pitch  angle  greater  than  90  degrees  for  values  of  6s 
less  than  about  -5  degrees.  This  is,  of  course,  true  up  to  the  pitchfork  bifurcation  point. 
This  phenomenon  of  inverted  pendulum  stabilization  was  first  discussed  by  McKinley 
[Ref.l4].  It  appears,  however,  from  the  results  of  Figure  18  that  the  inverted  pendulum 
stabilization  is  highly  degenerate  and  is  destroyed  as  soon  as  a  non-zero  rudder  angle  6r 
is  applied.  The  two  solutions,  6  and  ir-O,  do  not  cross  each  other  as  they  did  for  6r=0. 
Instead,  they  veer  of  each  other,  the  lower  one  corresponding  to  *=0  and  the  upper  one 
to  *=x.  As  a  result,  the  vehicle  configuration  is  always  keel  down,  the  only  difference 
between  the  two  solution  sets  at  Figures  17  and  18  being  180  degrees  heading. 
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Another  important  observation  can  be  made  by  considering  the  results  of  Figures 
16  and  17,  in  terms  of  the  roll  angle  «l».  It  appears  that  there  exists  a  range  of  stemplane 
angle  where  the  exact  value  of  $  is  very  sensitive  to  the  exact  value  of  5r  for  small 
values  of  6r.  The  rudder  angle  is  related  to  the  angle  of  attack  of  the  fluid  flow,  which 
is  highly  dependent  on  such  factors  as  vortex  shedding,  flow  separation,  and  ambient 
flow  turbulence  Saipkaya  [Ref,  5].  Since  these  factors  are  to  a  great  extent  uncertain  and 
difficult  to  model,  situations  like  the  one  described  above  should  be  avoided  in  practice. 
This  bifurcation  and  continuation  methods  used  in  this  thesis  provide  a  systematic  and 
effective  way  of  identifying  precisely  parameter  regions  where  such  sensitivity  of  the 
response  to  parameter  values  exists.  These  regions  should  be  avoided  when  executing  a 
proper  recovery  procedure. 
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rv.     CONTINUATION  IN  6r 

A.  GENERAL 

In  this  chapter,  the  results  of  continuation  in  terms  of  the  rudder  angle  will  be 
discussed.  In  this  continuation  run,  the  initial  conditions  were  taken  from  the  previous 
results  of  continuation  in  terms  of  the  rudder  angle.  The  configuration  of  the  vehicle  was 
just  like  the  previous  case.  In  this  case  the  system  of  equations  (35)  was  solved  for 
changing  values  of  the  primary  parameter  X  (here  X=  6r).  At  the  end  of  the  runs,  a 
hysteresis  behavior  was  observed  for  some  of  the  state  variables.  This  hysteresis 
behavior,  which  is  complementary  to  pitchfork  bifurcation,  is  explained  in  the  following 
section. 

B.  HYSTERESIS  PHENOMENON 

When  one-parameter  families  of  equation  are  studied,  there  is  no  need  to  talk  about 
a  hysteresis  phenomenon.  But  similarly  to  the  pitchfork  bifurcation,  this  situation  might 
change  when  a  two-parameter  family  of  equations  is  studied, 

f{y,X,a)     . 

Now  let  us  consider  the  foUowing  equation: 
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x^-X  +  ax  =  0 


(41) 


Let  us  call  this  algebraic  equation  as  a  function  of  x  F(x)  =  0.  It  can  be  observed  that 
equation  (41)  is  similar  to  the  generic  pitchfork  equation  (40).  The  only  difference  is  that 
the  roles  of  X  and  a  have  been  interchanged.  The  constant  term  is  now  our  primary 
bifurcation  parameter,  while  the  coefficient  of  x  is  the  perturbation  parameter.  We  know 
that  at  the  turning  point  F,  and  dF/dx,   should  be  zero.  So, 


df 
dx 


=0 


x^  =  -  — 


If  we  substitute  this  value  into  Equation  (41),  we  will  get  the  following  equation. 


27X2+4a3=0 


(42) 


This  equation  gives  the  relationship  between  X  and  a  at  the  vertical  tangent.  Now  we 


Figure  19.  Cusp  Curve 
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consider  equation  (41)  as  a  two  parameter  model  with  X  and  a  having  equal  rights. 
Equation  (42)  establishes  a  curve  in  the  (a  ,  X)  plane,  which  takes  the  form  of  a  cusp 
(Figure  19),  (Cusp  curves  will  be  explained  in  a  more  general  context  in  Chapter  VHI.), 
For  all  combinations  of  (  a,  X)  values  inside  the  cusp  (hatched  region),  Equation 
(41)  has  three  solutions;  for  values  outside  the  cusp,  there  is  only  one  solution.  This 
becomes  clear  when  branching  diagrams  are  drawn.  As  a  generic  behavior,  if  we  draw 
branching  diagrams  for  different  values  of  a  when  X  is  the  primary  continuation 
parameter,  we  get  the  Figures  20,  21,  22  according  to  the  sign  of  a. 

•  a  =  0 

•  a  >  0 

•  a  <  0 


X 


a  =  0 


Figure  20.    Branching  Diagram  When  o;  =  0 
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a  >  0 


Figure  21.    Branching  Diagram  When  a  <    0 


Figure  22 .    Branching  Diagram  When  a   >   0 

In  this  last  case,  we  see  a  hysteresis  effect.  When  hysteresis  occurs  in  a  branch, 
there  are  two  associated  turning  points  with  hysteresis  behavior.  These  points  occur 
precisely  at  the  two  real  solutions  of  X  from  Equation  42.  Here  there  is  a  jump  when  the 
solution  is  on  an  appropriate  part  of  the  branch.  The  two  outer  solution  in  x  are  stable 
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while  the  inner  most  solution  between  the  two  turning  points  is  unstable.  The  two  dotted 
straight  lines  of  Figure  19  represent  two  distinct  paths  through  the  cusp  and  explained 
in  Chapter  Vm. 

C.      STEADY  STATE  RESULTS 

Figures  23  through  30  show  steady  state  solutions  for  the  eight  state  variables, 
namely;  u,  v  ,w,  p,  q,  r,  «l>,  6.  In  Figures  23  through  30  solid  lines  correspond  to  5s  = 
0  degree,  dashed  lines  to  6s  =  -5  degrees,  dotted  lines  to  6s  =  -10  degrees  and  dash-dot 
lines  to  6s  =  -20  degrees.  For  these  continuation  runs,  the  rudder  angle  (6r)  is  the 
primary  continuation  parameter  (i.e,X)  and  stem  plane  angle  (6s)  is  perturbation 
parameter  (i.e,a)  and  values  of  all  other  parameters  are  kept  fixed:excess  buoyancy, 
6B  =  2  %  of  the  vehicle  weight  (W);  deflection  of  bow  planes,  6b  =  0  ;  location  of 
horizontal  and  vertical  centers  of  buoyancy  ,Xb=  Zg  =  0; location  of  longitudinal  center 
of  buoyancy,  x^b  =-1  %  of  the  vehicle  length  (L). 
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Figure  23.    Steady  State  Surge  Velocity 
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Figure  24.    Steady  State  Sway  Velocity 
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HEAVE  VELOCITY  W  vs  DR 
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-15 


-10 


-5 


RUDDER  ANGLE 

Steady  State  Heave  Velocity 


10 


1  5 


20 


ROLL  ANGULAR  VELOCITY  P  vs  OR 
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PITCH  ANGULAR  VELOCITY  U  vs  DR 
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Figure  27 .    Steady  State  Pitch  Angular  Velocity 
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Figure  28.    Steady  State  Yaw  Angular  Velocity 
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ROLL  ANGLE  PHI  vs  DR 
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Figures  23  through  30  show  that  the  horizontal  plane  state  variables  v,  p,  r,  $ 
behave  like  Figure  22  and  exhibit  a  hysteresis  behavior  for  certain  values  of  6s.  If  we 
zoom  in  Figure  28,  we  obtain  Figure  31  and  can  see  this  hysteresis  effect  very  clearly. 
When  6s  =  -15  degrees,  there  is  a  typical  hysteresis  phenomenon.  Figures  24,  26,  and 
28  also  show  the  same  kind  of  behavior. 
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Figure  31.    Typical  Hysteresis  Effect  from  Continuation 


On  the  other  hand,  the  vertical  plane  state  variables  u,  w,  q,  Q  behave  like  an  even 
function  in  x  (for  example  |  x  |  ).  This  behavior  can  be  seen  in  Figures  23,  25,  27,  30. 
where  it  is  clear  that  they  exhibit  generic  behavior  of  an  even  function  for  different 
values  of  the  perturbation  parameter  (i.e,  6s).  If  we  zoom  in  Figure  23,  we  can  see  this 
behavior,  very  clearly. 
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V.     CONTINUATION  IN  x 


GB 


A.  GENERAL 

In  this  chapter,  the  results  of  continuation  in  terms  of  the  longitudinal  center  of 
buoyancy  will  be  discussed.  In  continuation  runs,  the  system  of  equations  (35)  was 
solved  for  changing  values  of  the  primary  parameter  X  where  in  this  case  X  =  Xob).  In 
these  runs,  first  we  fixed  the  value  of  rudder  angle  and  changed  the  values  of  stem  plane 
angle,  then  we  kept  the  stemplane  angle  fixed  and  changed  the  rudder  angle,  to  see  how 
the  steady  state  values  were  changing.  When  we  keep  the  rudder  angle  at  zero  and 
change  the  values  of  stemplane  angle,  we  observe  that  there  is  a  pitchfork  bifurcation  for 
the  states  v,  p,  r.  Then,  we  fixed  8s  at  -20  degrees  and  started  to  change  mdder  angle 
to  see  how  the  pitchfork  bifurcation  diagrams  were  changing  with  changing  mdder  angle. 
As  a  last  continuation  mn  for  this  configuration,  the  value  of  stemplane  angle  was  fixed 
at  zero  and  we  changed  the  value  of  mdder  angle.  We  didn't  see  any  solution  branching 
for  this  case. 

B.  STEADY  STATE  RESULTS 

1.       5r=0  Degrees,  6s  is  Perturbation  Parameter 

Figures  33  through  40  show  steady  state  solutions  for  eight  state  variables, 
namely;  u,  v,  w,  p,  q,  r,  *,  6.  In  Figures  33  through  40,  solid  lines  correspond  to  8s=0 
degrees,  dotted  lines  correspond  to  5s =-5  degrees,  and  dash-dot  lines  correspond  to  6s=- 
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20  degrees.  In  graphs  of  the  state,  variables  v,  p,  r,  6  dash-dot  lines  may  seem  like  solid 
lines;  this  is  because  the  continuation  between  the  parameter  values  Xgb  =  -13  degrees 
and  Xgb  =  -7  degrees  traced  the  same  curve  twice  in  order  to  make  sure  that  all  solutions 
for  all  the  state  variables  were  captured.  Therefore,  out  of  plane  solutions  in  this 
parameter  range  belong  to  stemplane  angle  6s  =  -20  degrees. 

In  all  these  cases,  the  longitudinal  center  of  buoyancy  is  the  continuation  parameter 
and  6s  is  the  perturbation  parameter,  while  6r  is  kept  at  zero  degrees.  Values  of  the  other 
parameters  are  kept  fixed:  excess  buoyancy,  SB  =2  %  of  the  vehicle  weight  (W); 
deflection  of  the  bow  planes,  6b  =  0:  location  of  horizontal  and  vertical  centers  of 
gravity,  Xb  =  Zg  =0:  location  of  vertical  center  of  gravity,  Zgb=    1  feet. 
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ROLL  ANGULAR  VELOCITY  P  vs  XGB 
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Figure  36.    Steady  State  Roll  Angular  Velocity 
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Figure  37.    Steady  State  Pitch  Angular  Velocity 
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YAW  ANGULAR  VELOCITY  R  vs   XGB 
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Figure  38.    Steady  State  Yaw  Angular  Velocity 
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2.       5s=-20  Degrees,  5r  is  the  Perturbation  Parameter 

Figures  41  through  48  show  steady  state  solutions  for  the  state  variables  u, 
V,  w,  p,  q,  r,  4>,  6.  In  Figures  41  through  48,  solid  lines  correspond  to  6r=0  degrees, 
dotted  lines  correspond  to  6r=2.5  degrees,  and  dash-dot  lines  correspond  to  6r=5 
degrees. 

In  all  of  these  cases,  longitudinal  center  of  buoyancy  is  the  primary 
continuation  parameter  and  6r  is  the  perturbation  parameter,  while  6s  is  kq)t  at  -20 
degrees.  Values  of  the  other  parameters  are  kept  fixed:  Excess  buoyancy,  5B  =  2  %  of 
the  vehicle  weight  (W):  deflection  of  bow  planes,  6b  =  0;  location  of  horizontal  and 
vertical  centers  of  buoyancy,  Xb  =  Zg  =0;  location  of  the  vertical  center  of  gravity,  Zg= 
.1  feet. 
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Figure  43.    Steady  State  Heave  Velocity 


1.5 


ROLL  ANGULAR  VELOCITY  P  vs  XGB 


0.06 


0.04.- 


0.02  - 


Q.  0 


-0.02  - 


-0.04  - 


-0.06 


.'*'■''  ' .^^ 

\  •  •  ■ 

\ 
1 

■■■.\ 

V 

■    I                          ■  -*  ^..  - 

/ 
/■ 

«'•■•■' 

\ 

• 

'^L-.- •••.'■• 

,  ■ 

■  •j,;:-'=r:. 

^::::^'-- 

1 

-2 


Figure   44. 


-1.5 


-1 


-0.5 


0.5 


LONGITUDINAL  CENTER  OF  BUOYANCY  55L 

Steady  State  Roll  Angular  Velocity 


1.5 


56 


PITCH  ANGULAR  VELOCITY  Q  vs  XGB 


0.03 


0.025  - 


0.02  - 


O  0.015  - 


0.01  - 


0.005  - 


-2 


'.                          '■                          !,/■'■•. 

:              /'  :    .     /\.  /  :  ,        \ 

<.•  ■  •  •^ 

',  \ 
\             '.       !        \ 

\                                ^ 
\           ■    ■                  ^ 

/!.'/.■>                                              X 
1        '.         .'         1      .'    l'.                                                          V 
/        '.     .'        1    .'  1  '.                                                \ 
/            '.  ,'        1     .'  1     '.                                                  \ 

\  V  '■-., 

s 

.     1      '.   1            ■     y   .<              •                          ' 

:    X 

■••... 

'"*•*., 

-1.5 


-1 


-0.5 


0.5 


LONGITUDINAL  CENTER  OF  BUOYANCY  r,L 

Figure  45.    Steady  State  Yaw  Angular  Velocity 


1.5 


-2 


Figure  46 


YAW  ANGULAR  VELOCITY  R  vs   XGB 

0.04 

0.02 
0 

1 

•  /           ■  '■             \      • 

,.'•' 

■  ■  /■  ^■. :  ■  ■■:  ■  v;  •  •  "V : > 

1      \                     \.              / 

/          ^        ■         '  '           \                X 

/                \     !        /  -1              \              / 

V-..;-   ■.     :    \/ 

Y^-''^-::::: 

■'■'.::::::• 

,.v';\      '\;/\:>''' 

J 

/ 

0.02 

/ 
/ 

J         \ 

^^<^^:^ 

. .  •' 

iSs,^   ;  ^y^         .  *  s,     /        \ 

0.04 

V  ^^    '  ^^              '         ^>' 

V  ■  .     !               .  '         '.     / 

\     ■•!...-  ■                '..• 

\    '                         '■ 

s.                       ^      . 

n  OR 

-1.5 


-1 


-0.5 


0.5 


LONGITUDINAL  CENTER  OF  BUOYANCY  %L 

Steady  State  Yaw  Angular  Velocity 


1.5 


57 


YAW  ANGLE  PHI  vs   XGB 

50 
0 

'■                     / 

•                    /'■            ■■' 

:                        !   \              ■ 

\                     /     \ 

i /. ....ui. 

— 

•> 

— 

; 1 y.. 

.A '}      y 

•  '              ;,•'                   jy    : 
"v.    :        y'               .■•:                 yT 

-4^.  ..■•••■'  \  /     1 

•  ■              / 
I                                             / 

■                                 •        / 

•  :     / 

-2 


-1.5 


-1 


-0.5 


0.5 


1.5 


LONGITUDINAL  CENTER  OF  BUOYANCY  %L 
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3.       5s  =  0  Degrees,  5r  is  the  Perturbation  Parameter 

Figures  49  through  56  show  steady  state  solutions  for  the  eight  state 
variables,  namely  u,  v,  w,  p,  q,  r,  4>,  $.  In  Figures  49  through  56,  solid  lines 
correspond  to6r  =  0  degrees,  dashed  lines  correspond  to  6r  =  0.5  degrees,  dotted  lines 
correspond  to  6r  =  2.5  degrees  and  dash-dot  lines  correspond  to  6r  =  5  degrees. 

In  all  these  cases,  the  longitudinal  center  of  buoyancy  is  the  primary 
continuation  parameter  and  6t  is  the  perturbation  parameter.  Values  of  the  other 
parameters  are  kept  fixed:  Stemplane  action,  6s  =  0  degrees: excess  buoyancy,  5B  =  2 
%  of  the  vehicle  weight  (W);  deflection  of  bow  planes,  6b  =  0;  location  of  horizontal 
and  vertical  centers  of  buoyancy,  Xg  =  Zb  =  0;  location  of  vertical  center  buoyancy,  Zq^ 
=  0.1  feet. 
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Figure   50. 
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ROLL  ANGULAR  VELOCITY  P  vs  XGB 
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VI.     CONTINUATION  IN  5B 

A.      GENERAL 

In  this  chapter,  the  results  of  continuation  in  terms  of  the  excess  buoyancy  will  be 
discussed.  In  continuation  runs,  the  system  of  equations  (35)  was  solved  for  changing 
values  of  the  primary  parameter  X  (in  this  case  X=  SB).  In  these  runs,  first  we  fixed  the 
value  of  rudder  angle  and  changed  the  values  of  stemplane  angle,  then  we  kept  the 
stemplane  angle  fixed  at  zero  and  changed  the  rudder  angle  to  see  how  the  steady  state 
values  were  changing.  When  we  kept  the  rudder  angle  at  zero  and  changed  the  values  of 
stemplane  angle,  we  saw  that  there  was  a  pitchfork  bifurcation  for  the  state  variables  v, 
p,  r,  and  $.  After  we  observed  this  pitchfork  bifurcation  for  6s  =  -20  degrees  and  8r=0 
degrees,  we  started  to  change  the  rudder  angle  8r  to  see  how  the  pitchfork  bifurcation 
diagrams  were  changing.  As  an  other  continuation  run  we  fixed  6r=0  degrees  and 
changed  8s  towards  positive  values;  however,  we  didn't  see  any  branching  in  this  case. 
As  a  last  continuation  run  for  this  configuration,  the  value  of  stemplane  angle  was  fixed 
at  zero  and  the  value  of  mdder  angle  was  changed.  We  didn't  see  any  branching  behavior 
for  this  case. 
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B.       STEADY  STATE  DIAGRAMS 

1.       dr=0  Degrees,  ds  is  the  Perturbation  Parameter 

Figures  57  through  64  show  steady  state  solutions  for  eight  state  variables, 
namely;  u,  v,  w,  p,  q,  r,  4>,  6.  In  Figures  57  through  64,  solid  lines  correspond  to  6s =0 
degrees,  dashed  lines  correspond  to  6s =-5  degrees,  and  dash-dot  lines  correspond  to 
6s =-20  degrees.  In  all  of  these  cases,  excess  buoyancy  is  the  continuation  parameter, 
while  6r  is  kept  at  zero  degrees.  Values  of  the  other  parameters  are  kept  fixed: 
longitudinal  center  of  buoyancy  Xgb=-1  %  of  the  vehicles  length  (L);  deflection  of  the 
bow  planes,  6b =0;  location  of  horizontal  and  vertical  centers  of  buoyancy,  Xb=Zb=0; 
location  of  vertical  center  of  buoyancy,  Zgb=1  feet. 

As  we  can  see  from  Figures  57  through  64,  there  is  a  pitchfork  bifurcation  for  state 
variables  u,  v,  p,  r,  *.  We  have  out  of  plane  solutions  for  stemplane  angle  6s  =  -20 
degrees,  this  is  also  the  configuration  in  which  we  observe  bifurcations. 
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ROLL  ANGULAR  VELOCITY  P  vs  DELB 
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Figure  60.    Steady  State  Roll  Angular  Velocity 
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YAW  ANGULAR  VELOCITY  R  vs  DELB 
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After  analyzing  the  above  diagrams,  we  decided  to  look  at  a  case  with  a 
positive  sternplane  angle.  We  observed  that  there  was  no  solution  branching  for  this 
case.  We  got  the  inplane  solutions  and  there  were  no  out  of  plane  solutions.  In  Figure 
65,  we  can  see  the  situation  for  the  state  variable  u.  In  this  figure,  solid  lines  correspond 
to  6s=0  degrees,  dotted  lines  correspond  to  5s =5  degrees,  and  dash-dot  lines  correspond 
to  6s  =  20  degrees.  The  steady  state  solutions  are  not  very  sensitive  to  positive  sternplane 
action. 
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SURGE  VELOCITY  U  vs  DELB 
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Figure  65.    Steady  State  Surge  Velocity  When  6s  is  Positive 


2.       5s  =  -20  Degrees,    Br  is  the  Perturbation  Parameter 

Figures  66  through  73  show  steady  state  solutions  for  the  state  variables  u, 
V,  w,  p,  q,  r,  4»,  $.  In  these  figures,  dashed  lines  correspond  to  6r=  0.5  degrees,  and 
dash-dot  lines  correspond  to  6r=  5  degrees.  In  all  of  these  the  cases,  excess  buoyancy 
is  the  primary  continuation  parameter  and  6t  is  the  perturbation  parameter,  while  6s  is 
kept  -20  degrees.  Values  of  the  other  parameters  are  kept  fixed.  Longitudinal  center  of 
buoyancy,  -1  %  of  the  vehicle  length  (L);  deflection  of  bow  planes,  5b  =  0;  location  of 
horizontal  and  vertical  centers  of  buoyancy,  Xb=  Zb=  0;  location  of  the  vertical  center 
of  buoyancy,  Zgb=  0.1  feet.  By  doing  this,  we  saw  how  the  bifurcation  points  were 
changing  as  we  changed  the  perturbation  parameter  5r. 
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Figure  66.    Steady  State  Surge  Velocity 
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Figure  67.    Steady  State  Sway  Velocity 
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Figure  68.    Steady  State  Heave  Velocity 
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Figure  69.    Steady  State  Roll  Angular  Velocity 
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Figure  70.    Steady  State  Pitch  Angular  Velocity 
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Figure  71.    Steady  State  Yaw  Angular  Velocity 
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ROLL  ANGLE  PHI  vs  DELB 
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Figure  72.    Steady  State  Roll  Angle 
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Figure  73.    Steady  State  Pitch  Angle 
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3.       5s =0  Degrees,  dr  is  Perturbation  Parameter 

Figures  74  tlirough  81  show  steady  state  solutions  for  the  eight  state 
variables,  namely  u,  v,  w,  p,  q,  r,  4>,  6.  In  Figures  74  through  81  solid  lines  correspond 
to  5r=  0  degrees,  dashed  lines  correspond  to  6r=  0.5  degrees,  and  dash-dot  lines 
correspond  to  6r=5  degrees.  In  all  these  cases,  the  excess  buoyancy  SB  is  the  primary 
continuation  parameter  and  6r  is  the  perturbation  parameter.  Values  of  the  other 
parameters  are  kept  fixed:  Stemplane  action,  8s =0  degrees;  longitudinal  center  of 
buoyancy  \qq=-\%  of  the  vehicle  length  (L);  deflection  of  bowplanes,  5b =0;  location 
of  horizontal  and  vertical  centers  of  buoyancy,  Xb=  Zb=  0;  location  of  vertical  center  of 
buoyancy,  Zgb=  0.1  feet.  As  we  can  see,  there  was  no  solution  branching  for  this  case. 
Of  course,  we  got  out  of  plane  solutions  for  nonzero  rudder  angle  values. 
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Figure  74.    Steady  State  Surge  Velocity 
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Figure  76.    Steady  State  Heave  Velocity 
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Figure  77.    Steady  State  Pitch  Angular  Velocity 
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Figure  78.    Steady  State  Roll  Angle 
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Figure  79.    Steady  State  Yaw  Angular  Velocity 
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PITCH  ANGLE  THETA  vs  DELB 
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Figure  81.    Steady  State  Pitch  Angle 
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Vn.     CONTINUATION   IN  z 


GB 


A.  GENERAL 

In  this  chapter,  the  results  of  continuation  in  terms  of  the  vertical  center  of 
buoyancy  will  be  discussed.  In  continuation  runs,  the  system  of  equations  (35)  was  solved 
for  changing  values  of  the  primary  parameter  X  (in  this  case  X=  Zgb)-  In  these  runs,  first 
we  fixed  the  value  of  rudder  angle  and  changed  the  values  of  stemplane  angle,  then  we 
kept  the  stemplane  angle  fixed  at  zero  and  changed  the  rudder  angle  to  see  how  the 
steady  state  values  were  changing.  When  we  kept  the  rudder  angle  at  zero  and  changed 
the  values  of  stemplane  angle,  we  saw  that  there  was  a  pitchfork  bifurcation  for  the  state 
variables  u,  v,  w,  p,  r,  *,  6.  After  we  saw  that  there  was  a  pitchfork  bifurcation  when 
6s  =  -20  degrees  and  6r=  0  degrees,  we  started  to  change  mdder  angle  8t  to  see  how  the 
pitchfork  bifurcation  diagrams  were  changing.  As  a  last  continuation  mn  for  this 
configuration,  the  value  of  stem  plane  angle  was  fixed  at  zero  and  the  value  of  mdder 
angle  was  changed.  We  didn't  observe  any  branching  behavior  for  this  case. 

B.  STEADY  STATE  RESULTS 

1.       6r=  0  Degrees,  ds  is  the  Perturbation  Parameter 

Figures  82  through  89  show  steady  state  solutions  for  eight  state  variables, 
namely;  u,  v,  w,  p,  q,  r,  f>,  6.  In  Figures  82  through  89,  solid  lines  correspond  to  6s= 
0  degrees,  dashed  lines  correspond  to  6s =-5  degrees  and  dash-dot  lines  correspond  to 
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6s =-20  degrees.  In  all  of  these  cases,  the  vertical  center  of  buoyancy  is  the  primary 
continuation  parameter,  and  5s  is  the  perturbation  parameter,  while  6r  is  kept  at  zero 
degrees.  Values  of  the  other  parameters  are  kept  fixed:  Excess  buoyancy  5B=2%  of  the 
vehicle  weight  (W);  longitudinal  center  of  buoyancy  \q^=-1%  of  vehicle  length; 
deflection  of  the  bow  planes  8b =0;  location  of  the  vertical  and  horizontal  centers  of 
buoyancy,  Zb=Xb=  0.  As  we  see  from  Figures  82  through  89,  there  is  a  pitchfork 
bifurcation  for  state  variables  u,  v,  w,  p,  r,  4»,  6.  We  have  out  of  plane  solutions  for 
this  vehicle  configuration. 
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SWAY  VELOCITY  V  vs  ZGB 
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ROLL  ANGULAR  VELOCITY  P  vs  2GB 
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YAW  ANGULAR  VELOCITY  R  vs  ZGB 
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2.       5s=-20  D^rees,  hv  is  the  Perturbation  Parameter 

Figures  90  through  97  show  steady  state  solutions  for  the  state  variables  u, 
V,  w,  p,  q,  r,  4»,  B.  In  Figures  90  through  97,  solid  lines  correspond  to  6r=0  degrees, 
dashed  lines  when  6r=0.5  degrees,  dash-dot  lines  correspond  to  6r=5  degrees.  In  all 
of  these  cases,  vertical  center  of  buoyancy  Zqb  is  primary  continuation  parameter  and  6r 
is  the  perturbation  parameter,  while  6s  is  kept  at  -20  degrees.  Values  of  the  other 
parameters  are  kept  fixed;  excess  buoyancy  8b =2%  vehicle  weight  (W),  longitudinal 
center  of  buoyancy  Xgb--1  %  of  the  vehicle  length  (L);  deflection  of  bowplanes  5b =0; 
location  of  horizontal  and  vertical  centers  of  buoyancy,  Xb=Zb=0.  By  doing  this,  we 
saw  how  the  bifurcation  points  were  changing  as  we  change  perturbation  parameter  8r. 
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A  bifurcation  phenomenon  took  place  for  6r=0  degrees  and  disappeared  after  we  applied 
a  small  nonzero  rudder  angle  value. 
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ROLL  ANGLE  PHI  vs  ZGB 
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3.       5s  =  0  Degrees,  Br  is  the  Perturbation  Parameter 

Figures  98  through  105  show  steady  state  solutions  for  the  eight  state 
variables,  namely  u,  v,  w,  p,  q,  r,  *,  6.  In  Figures  98  through  105,  solid  lines 
correspond  to  6r=0  degrees,  dashed  lines  correspond  to  6r=.5  degrees,  dash-dot  lines 
correspond  to  8r=5  degrees.  In  all  these  cases,  the  vertical  center  of  buoyancy  is  the 
primary  continuation  parameter.  Values  of  the  other  parameters  are  kept  fixed;  stemplane 
action  8s =0  degrees;  longitudinal  center  of  buoyancy  x^b  =  -1  %  of  vehicle  length  (L); 
deflection  of  mow  planes,  8b =0  degrees;  excess  buoyancy  8B=  2%  of  vehicle  weight 
(W);  location  of  horizontal  and  vertical  centers  of  buoyancy  Xb=Zb=0.  As  we  see,  there 
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is  no  solution  branching  for  this  last  case.   We  have  out  of  plane  solutions  for  nonzero 
rudder  angles. 
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Vm.     CUSP  CURVES 

A.      GENERAL 

In  this  chapter,  the  cusp  phenomenon  will  be  discussed.  After  seeing  the  branching 
characteristics  of  the  steady  state  solutions  of  our  system,  we  decided  to  look  at  how  the 
bifurcation  points  and  the  turning  points  were  changing  with  changing  values  of  rudder 
angle.  The  main  difference  between  a  bifurcation  point  and  a  turning  point  is  that,  at  a 
bifurcation  point,  exactly  two  branches  intersect  with  two  distinct  tangents;  whereas,  at 
a  turning  point,  the  branch  comes  from  one  side  and  turns  back.  We  can  see  this 
phenomenon  from  Figure  106.  This  is  the  graph  that  was  generated  for  the  surge  velocity 
when  the  primary  continuation  parameter  was  x^Band  the  stemplane  angle  6s =-20 
degrees.  In  Figure  106,  points  A  and  B  are  bifurcation  points,  and  points  C  and  D  are 
turning  points.  Clearly,  locally  there  are  no  solutions  on  one  side  of  a  turning  point  and 
two  solutions  on  the  other  side,  but  there  are  solutions  on  each  side  of  a  bifurcation 
point. 

Let  us  look  at  Equation  43  again,  this  is  the  model  equation  that  represents  a  cusp, 

x^-Bx+A=Q  (43) 

Let  us  call  this  algebraic  equation  as  a  function  of  x,  F(x)=0.  For  the  critical  curve,  F 
and  dF/dx  should  be  zero. 
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-^=0       =>       x^  =  l 
dx  3 


If  we  substitute  this  value  into  equation  (43),  we  get  the  following  relation: 

If  we  plot  this  relation  in  the  AB  plane,  we  get  the  cusp  curve  shown  in  Figure  107.  This 
represents  a  separation.  Equation  43  has  three  solutions  for  parameter  values  inside  the 
cusp  and  one  solution  outside  this  parameter  values.  In  Figure  107,  if  we  fix  the  value 
of  A  and  vary  the  value  of  B  (path  (a)  in  Figure  107)  we  get  a  pitchfork  bifurcation 
diagram,  as  shown  in  Figure  108. 

CONSTANT  A,   CHANGING   B 


Figure  108.    Pitchfork  Bifurcation 

In  Figure  107,  if  we  fix  the  value  of  the  value  of  B  and  vary  the  value  of  A  (path  (b)  in 

Figure  107)  we  get  a  hysteresis  behavior,  as  shown  in  Figure  110.  It  is  clear,  that  when 
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a  line  crosses  the  cusp  curve  in  the  (A,  B)  plane  a  branching  phenomenon  takes  plane. 
As  mentioned  before,  Equation  43  has  three  solutions  for  parameter  values  inside  the 
cusp  and  one  solution  outside  the  cusp  this  can  be  seen  easily  in  Figure  109.  Now  let  us 
take  a  look  at  Figure  111.  As  can  be  seen,  there  is  a  simple  way  of  recovering  a 
bifurcation  diagram  from  the  path  through  the  given  cusp  by  lifting  the  path  to  the  cusp 
surface. 


Figure  109.    3-D  Graph  of  Equation  43 


99 


CONSTANT  B,  CHANGING  A 


-1  -0.8        -0.6        -0.4        -0.2 


0.2  0.4.  0.6  0.8  1 


Figure  110.    Hysteresis  Diagram 


Figure  111.    Geometry  of  Cusp  Phenomenon 
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If  we  have  the  cusp  curves  for  a  multiparameter  equations  system,  we  can  always 
have  an  idea  about  the  solutions  of  this  system  for  different  perturbed  cases.  This  means 
that  by  looking  at  the  path  of  the  parameters  through  the  cusp  curve,  we  can  see  how 
many  solutions  we  have  for  the  corresponding  configuration  of  the  system.  Similarly,  we 
can  think  about  another  problem  which  has  a  turning  point.  If  we  draw  the  position  of 
turning  point  in  the  plane  with  changing  parameter  values,  we  will  get  two  solutions 
inside  the  curve  and  no  solutions  outside  these  parameter  values.  Therefore,  these  graphs 
can  tell  us  the  nature  of  the  solution  set  of  our  system.  In  our  study,  we  tried  to  get  the 
cusp  curves  and  the  position  of  the  turning  points  and  the  results  are  presented  in  the 
following  section. 

B.      RESULTS 

First  of  all,  let  us  look  at  the  case  when  we  have  the  stemplane  action  as  our 
primary  continuation  parameter.  When  6s  is  between  -5  and  -10  degrees,  we  had  a 
bifurcation  point  and  in  that  case  the  rudder  angle  was  zero,  see  Figure  13.  To  get  the 
complete  cusp  curve,  we  changed  the  value  of  6r  and  found  the  locations  of  the 
bifurcation  point.  These  were  plotted  as  6s  versus  6r,  see  Figure  112. 
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CUSP  CURVE  WHEN  CONTINUATION  PARAMETER  is  DS 
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Figure  112.     Cusp  Curve  When  Continuation  Parameter  is 
Sternplane  angle         _. 


In  Chapter  V,  we  saw  that  when  we  have  Xqb  as  primary  continuation  parameter, 
and  sternplane  action  5s =-20  degrees,  we  have  both  bifurcation  and  turning  points. 
Figure  106.  To  see  the  effect  of  changing  the  rudder  angle  on  the  bifurcation  and  turning 
points,  we  changed  the  value  of  rudder  angle  and  found  the  location  of  bifurcation  and 
turning  points  in  the  (Xqb,  dt)  plane.  The  results  are  shown  on  the  following  graphs.  The 
bifurcation  points  generated  a  cusp  curve  and  the  turning  points  generated  a  parabola  in 
the  (xBo,5r)  plane. 
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Figure  113.    Cusp  Curve  When  Continuation  Parameter  is  x^g 
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Figure  115.   Second  Cusp  Curve  VOien  Continuation  Parameter  is 
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In  Chapter  VI,  we  used  the  excess  buoyancy  as  our  primary  continuation 
parameter.  In  that  configuration,  we  had  a  turning  point  and  a  bifurcation  point  for 
stemplane  angle  5s =-20  degrees,  (Figure  57  in  Chapter  VI).  To  see  the  effect  of 
different  rudder  angles  on  these  points,  we  changed  the  value  of  the  rudder  angle  and 
found  the  location  of  these  points  in  the  (Xob,  5r)  plane.  These  results  are  shown  in 
Figure  117  and  Figure  118.  Again,  we  got  a  cusp  curve  and  a  parabola. 
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Figure  118.    Turning  Point  When  Continuation  Parameter  is 
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In  Chapter  Vn,  we  used  the  vertical  center  of  buoyancy  as  our  primary 
continuation  parameter.  In  that  continuation,  run  we  found  that  there  were  two 
bifurcation  points  at  stemplane  angle  6s =-20  degrees,  see  Figure  82  Chapter  VII.  In  that 
figure,  the  bifurcation  point  when  Zqb  is  around  zero  is  not  clear,  but  it  becomes  clear 
when  we  changed  the  value  of  rudder  angle.  To  construct  the  cusp  curves,  we  changed 
the  value  of  rudder  angle  and  found  the  coordinates  of  bifurcation  points  in  the  (Zqb,  6r) 
plane.  These  graphs  are  shown  in  Figures  119  and  120. 
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Figure  119.    Cusp  Curve  When  Continuation  Parameter  is  z 
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As  a  final  study  in  this  chapter,  we  changed  the  value  of  yo  in  the  first  continuation 
run  and  tried  to  get  the  biased  cusp  curves  for  this  perturbed  case.  As  a  result,  we  can 
still  see  a  symmetric  bifurcation  for  this  perturbed  case  for  a  specific  parameter 
combination.  It  is  clear,  however,  that  this  situation  is  very  degenerate  and  can  easily  be 
destroyed  for  different  parameter  values.  In  Figure  121,  solid  lines  correspond  to  yG=0 
feet,  dashed  lines  correspond  to  yc =0.001  feet,  dotted  lines  correspond  to  yc =0.003 
feet,  and  dash-dot  lines  correspond  to  yc =0.005  feet. 
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Figure  121.   Cusp  Curve  When  Continuation  Parameter  is  6s  for 
Different  Values  of  y^ 
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K.     DYNAMIC  STABILITY  ANALYSIS 

A.  GENERAL 

So  far,  we  have  performed  continuation  runs  and  have  obtained  steady  state 
solutions  of  the  system.  For  stability  analysis,  we  chose  a  representative  case  and 
performed  the  stability  analysis  for  that  case.  We  selected  the  case  when  we  had  x^b  as 
our  primary  continuation  parameter  and  rudder  angle  6r=0  degrees,  and  stemplane  angle 
6s =-20  degrees.  To  predict  dynamic  stability,  we  used  the  six  degree  of  freedom 
equations  of  motion  along  with  the  Euler  angle  rate  equations  for  the  derivatives  of  the 
angles  of  pitch  and  roll  (6  and  4).  These  equations  of  motion  were  then  linearized  around 
the  steady  state  nominal  points  computed  in  the  previous  chapters.  Eigenvalue  analysis 
was  then  used  to  compute  the  stability  characteristics  of  each  solution. 

B.  DYNAMIC  STABILITY  RESULTS 

The  stability  analysis  results  are  shown  in  Figure  122  and  Figure  123.  In  these 
figures,  solid  lines  correspond  to  stable  solutions  while  dashed  lines  correspond  to 
unstable  solutions.  In  Figure  122  up  to  point  A,  we  have  a  stable  solution,  which  is 
predominantly  forward  motion,  because  we  have  a  high  surge  velocity  and  the  vehicle 
is  moving  in  the  vertical  plane.  Between  point  A  and  B  an  out  of  plane  solution  exists 
and  the  vehicle  exhibits  motion  in  the  horizontal  plane,  also.  From  Figure  123,  it  is  clear 
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that  at  point  A  the  sway  velocity  becomes  nonzero.  This  out  of  plane  solution  is  stable 
up  to  point  B,  where 


10 


6  - 


2  - 


0  :: 


SURGE  VELOCITV 


■  ■— ^ 

'^::^-'^^^^^'^'''y^-'- 

■jj,::---^ 

\^i:z:iz 

.  S ./ , - 

1            .    \ 

Si.fc.^fc*—— • 

C 

-1.5  -1  -O.S  0  O.S 

XCB 


Figure  122.    Steady  State  Surge  Velocity 
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the  out  of  plane  motion  becomes  unstable  up  to  point  C,  where  it  disappears.  In  Figure 
122,  the  motion  of  the  vehicle  from  point  D  to  E  is  stable  and  has  a  low  surge  velocity, 
which  means  that  the  motion  of  the  vehicle  is  predominantly  in  the  upward  direction. 
The  vehicle  is  moving  in  the  vertical  plane  and  there  is  no  horizontal  plane  motion. 

To  assess  stability  of  each  solution,  we  used  eigenvalue  analysis.  The  eigenvalues 
that  we  computed  for  the  following  cases  are  summarized  in  Table  1 . 

•  Case  1 

Stemplane  angle  8s=-20  degrees,  rudder  angle  8r=0  degrees,  Xgb=-2  %  L,  in 
plane  solution. 

•  Case  2 

Stemplane  angle  6s =-20  degrees,  rudder  angle  6r=0  degrees,  Xgb=-1  25  %  L,  out 
of  plane  solution. 

•  Case  3 

Stemplane  angle  8s =-20  degrees,  rudder  angle  8r=0  degrees,  x<3b=0.57%  L  point 
C   Figure  122. 

•  Case  4 

Stemplane  angle  8s =-20  degrees,  mdder  angle  8r=0  degrees,  x<3b=0  %  L.  This 
is  on  the  stable  branch  between  point  D  and  E  in  Figure  122. 

Table  1  represents  the  eigenvalues  of  the  system.  We  can  see  for  case  three  that 
there  is  a  pair  of  complex  conjugate  eigenvalues  with  positive  real  parts.  This  means  that 
we  have  a  hopf  bifurcation  prior  to  that  point,  and  we  expect  periodic  solutions  for  that 
part  of  the  solution  branch.  Since  we  have  a  stable  solution  up  to  point  E  in  Figure  122, 
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we  don't  expect  to  see  this  periodic  out  of  plane  behavior  up  to  point  E.  The  methods 
that  we  used  in  this  work  give  only  stationary  solutions,  therefore,  to  see  this  periodic 
behavior  we  used  numerical  simulation  techniques. 

TABLE  1 


Eigenvalues 

Case  1 

Case  2 

Case  3 

Case  4 

X, 

-0.962+0.499 

-2.407 

-2.591 

-0.093+0.522 

K 

-0.962-0.499 

-0.855-h0.413 

-0.968 -»- 0.427 

-0.093-0.522 

X, 

-0.068-1-0.014 

-0.855-0.413 

-0.968-0.427 

-0.067+0.064 

X4 

-0.068-0.014 

-0.723 

-1.048 

-0.067-0.064 

X5 

-0.761 

-0.066 

-0.314 

-0.193+1.531 

K 

-2.731 

-0.018-1-0.044 

-0.129 

-0.193-1.531 

X7 

-0.198-1-0.115 

-0.181-0.044 

-f- 0.015 +0.004 

-0.240 

X, 

-0.198-0.115 

-0.271 

+0.015-0.004 

-0.045 

C.     NUMERICAL  INTEGRATION  RESULTS 

The  linearized  dynamic  response  results  were  verified  by  simulations  using 
numerical  integration  of  the  full  six  degrees  of  freedom  equations  of  motion  for  the 
swimmer  delivery  vehicle  (SDV).  Figures  124  and  125  show  a  plot  of  surge  velocity 
(U),  and  sway  velocity  (V)  versus  time  respectively  for  the  center  of  gravity  aft  of  the 
center  of  buoyancy  case  (Xob=-2)  with  a  diveplane  angle  (6s)  of  -20  degrees,  and  rudder 
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Figure  124.    Time  History  of  Surge  Velocity  for  Case  1 
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Figure  125.    Time  History  of  Sway  Velocity  for  Case  1 
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angle  (6r)  of  0  degrees.  The  steady  state  results  of  the  numerical  integration  method 
match  the  linearized  dynamic  results  exactly. 

Figures  126  and  127  show  a  plot  of  surge  velocity  (U),  and  sway  velocity  (V) 
versus  time  respectively  for  the  center  of  gravity  aft  of  the  center  of  buoyancy  case 
(XoB=  -1.25)  with  a  diveplane  angle  (6s)  of  -20  degrees  and  rudder  angle  (6r)  of  0 
degrees.  The  steady  state  results  of  the  numerical  integration  method  match  the  linearized 
dynamic  results  exactly.  Vehicle  has  an  out  of  plane  motion  in  this  case.  As  can  be  seen 
in  Figure  126.  the  sway  velocity  V  converged  to  a  nonzero  steady  state  value. 
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Figure  126.    Time  History  of  Surge  Velocity  for  Case  2. 
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Figure  127.    Time  History  of  Sway  Velocity  for  Case  2. 


Figures  128  and  129  show  a  plot  of  surge  velocity  (U),  and  sway  velocity  (V) 
versus  time  respectively  for  the  Xob=-0.75  %  L  with  a  diveplane  angle  (6s)  of  -20 
degrees,  and  rudder  angle  (6r)  of  0  degrees.  Here  we  can  see  the  periodic  out  of  plane 
behavior  of  the  steady  state  solution.  This  is  what  we  expected  from  the  eigenvalue 
analysis  that  we  did  in  the  previous  section. 

Figures  130  and  131  show  a  plot  of  surge  velocity  (U),  and  sway  velocity  (V) 
versus  time  respectively  for  Xgb=0  %  L  with  a  dive  plane  angle  (6s)  of  -20  degrees,  and 
rudder  angle  (6r)  of  0  degrees.  And  once  again,  the  steady  state  results  of  the  numerical 
integration  method  match  the  linearized  dynamic  results  exactly. 
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Figure  128.    Time  History  of  Surge  Velocity  for  Case  3. 
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Figure  129.    Time  History  of  Sway  Velocity  for  Case  3. 
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Figure  130.    Time  History  of  Surge  Velocity  for  Case  4 
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Figure  131.    Time  History  of  Sway  Velocity  for  Case  4 
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X.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  problem  of  steady  state  response  and  dynamic  stability  analysis  of  submarines 
in  free  positive  buoyancy  ascent  has  been  studied.  The  main  conclusions  of  this  work  can 
be  summarized  as  follows: 

•  Steady  state  motion  is  not,  in  general,  restricted  to  the  vertical  plane,  there  can 
exist  complicated  out-of-plane  motions.  A  combination  of  analytical  results  for  the 
in  plane  motions  and  numerical  simulations  and  continuation  was  used  in  order  to 
identify  all  possible  steady  states. 

•  For  some  steady  state  solutions  there  can  be  pitchfork  and  hysteresis  bifurcations. 
In  fact,  it  was  shown  that  the  pitchfork  is  the  primary  mechanism  for  generating 
out  of  plane  motions.  --^    — 

•  The  response  of  the  vehicle  can  be  very  sensitive  to  system  parameters.  In  such  a 
case,  any  numerical  integration  results  should  be  viewed  with  extreme  caution. 

•  The  biftircation  graphs  that  we  generated  are  very  useful  in  order  to  understand  the 
steady  state  solution  set  of  the  vehicle  for  any  parameter  combination. 

•  As  a  recommendation  for  future  research,  continuation  methods  should  be  used  to 
investigate  the  nature  of  the  Hopf  bifurcation,  specifically  stability  of  periodic 
solutions. 

•  Furthermore,  the  effect  of  nonzero  propulsion  should  be  investigated.  It  is  possible 
that  small  nonzero  propulsion  forces  could  affect  significantly  the  bifurcations. 
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APPENDIX 

1.       BIFPACK 

Biijpack  is  a  program  package  for  calculating  bifurcations  written  in 
FORTRAN  by  R.  Seydel  [Ref.  15] 

In  this  work  we  used  Hijack  to  find  the  bifurcation  characteristics  of  steady 
state  solutions  of  submersibles.  The  resulting  set  of  algebraic  equations  was  solved  by 
using  PACKA  which  is  devoted  to  system  of  nonlinear  equations 

f{y,k)=o 

Here  y  and  f(y,X)  are  vectors  with  n  components. 

To  use  Bifjpack  we  did  not  bother  about  how  to  program  our  original  f(y,X) 
into  the  enlarged  form  that  is  used  by  Bifipack.  This  larger  frame  that  embraces  the 
original  n  scalar  algebraic  equations  by  Bifipack  is  called  FRAMEA  in  Biiipack.  To  solve 
our  problem  we  duplicated  this  frame  and  gave  it  different  a  name.  This  became  our 

main  program  and  then  we  entered  n  scalar  formulas  f,,  f2  ,f3, ,  fn  into  subroutine 

FCN.  A  sample  of  those  two  programs  is  included  in  the  following  sections.  In 
conjunction  with  subroutine  FCN,  we  also  wrote  functions  which  compute  the  integral 
values  in  Equations  8,  9,  11,  and  12.  These  are  functions  FI2,  FI3,  FI5,  FI6  in 
subroutine  FCN. 
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2.       MAIN  PROGRAM 

C  SUB2.F0R 

PROGRAM  MAIN 
IMPUCIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  EPS 

EXTERNAL  FCN 

C 

C       INTERACTIVE  WORK  WITH  THE   PACK-A   SUBPACKAGE  OF 

C       BIFPACK 

C 

COMMON/DATA/DS,DB,DElJB,XGB,ZGB,XB,ZB,DRS,DRB 

OPEN  (11,  FILE  =  'DATAA1',STATUS  =  '0LD') 

OPEN  (8  ,  FILE  =  'DIAGRAM',STATUS  =  'OLD') 

OPEN  (3  ,  FILE  =  'START',STATUS  =  'OLD') 

OPEN  (1  ,  FILE  =  'OPTIONS', STATUS  = 'OLD') 

OPEN  (10,  FILE  ='FORTRAN.MAT',STATUS  =  'OLD') 

OPEN  (20,  FILE  ='FORTRANl. MAT', STATUS  = 'OLD') 

OPEN  (15,  FILE='DEN.MAT',STATUS  =  'OLD') 

OPEN  (16,  FILE='DEN1. MAT, STATUS  = 'OLD') 
C 

C   PROBLEM  DEPENDENT  VARIABLES: 

C   EPS  =  (REL.)  ACCURACY,   N  =  NUMBER  OF  EQUATIONS 
C 

EPS  =  l.E-5 

N=9 

JACOBI  =  l 
C 

C       DATA 
C 

PI      =4.0*ATAN(1.0) 

WEIGHT=  12000.0 

L=  17.425 

DB=0.0*PI/180.0 

DELB=2.0*WEIGHT/100.0 

XGB=-1.0*L/100.0 

ZGB=0.1 

XB=0.0*L/100.0 

ZB=0.0 

DRS=0.0*PI/ 180.0 
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DRB=0.0*PI/180.0 
WRITE  (*,*)    'WHAT  IS  DRS' 
READ  (*,*)  DRS 
DRS=DRS*PI/180.0 
WRITE  (*,*)   DRS 

WRITE  (11,*)  'DS  AS  PARAMETER  ZGB  =  .1,DB  IS  2  PERCENT 

OF  WEIGHT' 
WRITE(11,*)  '  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  -  = 


c 

c 

nSHTIAL  SOLUTION  USED  FO 

c 

c 

Y(l)  =  8.6 

U 

c 

Y(2)=0 

V 

c 

Y(3)=-0.25 

W 

c 

Y(4)=0 

P 

c 

Y(5)=0 

Q 

c 

Y(6)=0 

R 

c 

Y(7)=0 

PHI 

c 

Y(8)=0.53 

THETA 

c 

Y(9)=0 

PSIDOT 

WRITE  (*,*)'  Previous  data  in  DIAGRAM  to  be  purged?' 

WRITE  (*,*)  '      ENTER  0   FOR  PURGING,  ENTER    1    FOR 
APPENDING:' 

READ  (*,*)  NPURGE 

IF  (NPURGE.EQ.O)  REWIND  8 

IF  (NPURGE.EQ.1)  THEN 
DO  23  1=1,9999 
23  READ  (8,'(  )',END=29) 

29  CONTINUE 

END  IF 

CALL  INTERA  (N,FCN,EPS,JACOBI) 

WRITE  (*,*)  '  PRINT-OUTPUT  IN  FILE  "DATAAl"' 

STOP 

END 

3.       SUBROUTINE  FCN. 

C  EOMl.FOR 

SUBROUTINE  FCN  (TDUMMY,Y,F) 
IMPUCrr  DOUBLE  PRECISION  (A-H,0-Z) 
DOUBLE  PRECISION  Y,F,PAR,A,ETACO,ZETACO,PARCO,TDUMMY 


121 


DOUBLE  PRECISION  DB,  DR,  XG,  ZG,  YG,  XPROP,  12,  13,  YB, 
1  KPROP,  15,  XB,  ZB,  BOY,  16  ,NPROP,  U,  V,  W, 

2  FI2,  n3,n5,n6 

DOUBLE  PRECISION  WEIGHT,  IX,  lY,  IZ,  IXY,  lYZ,  IXZ,  L,  RHO, 

1  G,  AO,  CDO,  M,  XPP,  XQQ,  XRR,  XPR,  XUDOT, 

2  XWQ,  XVP,  XVR,  XQDS,  XQDB,  XDRDR,  XRDR,  XVV, 

3  XWW,XVDR,XWDS,XWDB,XDSDS,XDBDB,XRES,  YPDOT, 

4  YRDOT,YPQ,YQR,YVDOT,YP,YR,YVQ,YWD,YWR,YV, 

5  YVW,  YDR,  YDRB,  ZW,  ZDS,  ZDB,  KPDOT,  KVW, 

6  MQDOT,  MPP,  MPR,  ZQDOT,  ZPP,  ZPR,  ZRR, 

7  ZWDOT,ZQ,ZVP,ZVR,ZW,NRDOT,NPQ,NQR,  KRDOT, 

8  KPQ,  KQR,  KVDOT,  KP,  KR,  KVQ,  KWP,  KWR,  KV, 

9  MRR,  MWDOT,  MQ,  MVP,  MVR,  MW,  MVV,  MDS,  MDB, 

1  NPDOT,NVDOT,  NP,  NR,  NVQ,  NWP,  NWR,  NV,  NVW, 

2  NDR,  NDRB,  XRDRB,  XRDRS 
PARAMETER  (NDIM  =  12,NDIML=13) 
DIMENSION  Y(NDIML),F(NDIML) 

C 

C       FOLLOWING  THREE  STATEMENTS  ARE  FOR  CONTINUATION  AN 

C      BRANCH  SWITCHING  PURPOSES: 

C 

COMMON/DATA/DS,DB,DELB,XGB,ZGB,XB,ZB,DRS,DRB 

COMMON  /CONTRO/  ETACO,ZETACO,PARCO,NCO,NlCO,N2CO, 
1  KCO,IND2CO,NFLAG 

PAR=Y(NlCO) 

F(NlCO) = Y(KCO)-ZETACO*Y(IND2CO)-ETACO 
C 

C        GEOMETRIC  PROPERTIES 
C 

WEIGHT=  12000.0 

IX  =  1760.0 

lY =9450.0 

IZ= 10700.0 

IXY=0.0 

IYZ=0.0 

IXZ=0.0 

L=  17.425 

RHO  =  1.94 

G=32.2 

YG=0.1 

YB=0.0 

A0=2.0 

CD0=0.0057 
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MASS=WHGHT/G 

C 

C       SURGE  HYDRODYNAMIC  COEFHCIENTS 

C 

XPP     =  7.030E-03*0.5*RHO*L**4 
XQQ     =-1.470E-02*0.5*RHO*L**4 
XRR     =  4.010E-03*0.5*RHO*L**4 
XPR     =  7.640E-04*0.5*RHO*L**4 
XUDOT  =-7.580E-03*0.5*RHO*L**3 
XWQ     =-1.920E-01*0.5*RHO*L**3 
XVP     =-3.240E-03*0.5*RHO*L**3 
XVR     =  1.890E-02*0.5*RHO*L**3 
XQDS    =  2.610E-02*0.5*RHO*L**3 
XQDB    =-2.600E-03*0.5*RHO*L**3 
XRDR    =-8.180E-04*0.5*RHO*L**3 
XW     =  5.290E-02*0.5*RHO*L**2 
XWW     =  1.710E-01*0.5*RHO*L**2 
XVDR   =  1.730E-03*0.5*RHO*L**2 
XWDS    =  4.600E-02*0.5*RHO*L**2 
XWDB    =  9.660E-03*0.5*RHO*L**2 
XDSDS  =-1.160E-02*0.5*RHO*L**2 
XDBDB  =-8.070E-03*0.5*RHO*L**2 
XDRDR  =-1.010E-02*0.5*RHO*L**2 
XRES    =  CD0*0.5*RHO*L**2 

XPROP  =  XRES*ALPHA**2 
XRDRB=XRDR 
XRDRS=XRDR 
XVDRS=XVDR 
XDRB=XVDR 

C 

C       SWAY  HYDRODYNAMIC  COEFFICIENTS 

C 

YPDOT  =  1.270E-04*0.5*RHO*L**4 
YRDOT  =  1.240E-03*0.5*RHO*L**4 
YPQ     =  4.125E-03*0.5*RHO*L**4 
YQR     =-6.510E-03*0.5*RHO*L**4 
YVDOT  =-5.550E-02*0.5*RHO*L**3 
YP      =  3.055E-03*0.5*RHO*L**3 
YR      =  2.970E-02*0.5*RHO*L**3 
YVQ     =  2.360E-02*0.5*RHO*L**3 
YWP     =  2.350E-01*0.5*RHO*L**3 
YWR     =-1.880E-02*0.5*RHO*L**3 
YV      =-9.310E-02*0.5*RHO*L**2 
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YVW     =  6.840E-02*0.5*RHO*L**2 
YDRS    =+2.270E-02*0.5*RHO*L**2 
YDRB    =+2.270E-02*0.5*RHO*L**2 

C 

C       HEAVE  HYDRODYNAMIC  COEFFICIENTS 

C 

ZQDOT  =-6.810E-03*0.5*RHO*L**4 
ZPP     =  1.270E-04*0.5*RHO*L**4 
ZPR     =  6.670E-03*0.5*RHO*L**4 
ZRR     =-7.350E-03*0.5*RHO*L**4 
ZWDOT  =-2.430E-01*0.5*RHO*L**3 
ZQ      =-1.350E-01*0.5*RHO*L**3 
ZVP     =-4.810E-02*0.5*RHO*L**3 
ZVR     =  4.550E-02*0.5*RHO*L**3 
ZW      =-3.020E-01*0.5*RHO*L**2 
ZW     =-6.840E-02*0.5*RHO*L**2 
ZDS     =-2.270E-02*0.5*RHO*L**2 
ZDB     =-2.270E-02*0.5*RHO*L**2 

C 

C       ROLL  HYDRODYNAMIC  COEFHCIENTS 

C 

KPDOT  =-1.010E-03*0.5*RHO*L**5 
KRDOT  =-3.370E-05*0.5*RHO*L**5 
KPQ     — 6.930E-05*0.5*RHO*L**5 
KQR     =  1.680E-02*0.5*RHO*L**5 
KVDOT  =  1.270E-04*0.5*RHO*L**4 
KP      =-1.100E-02*0.5*RHO*L**4 
KR      =-8.410E-04*0.5*RHO*L**4 
KVQ     =-5.115E-03*0.5*RHO*L**4 
KWP     =-1 .270E-04*0.5*RHO*L**4 
KWR     =  L390E-02*0.5*RHO*L**4 
KV      =  3.055E-03*0.5*RHO*L**3 
KVW     =-L870E-01*0.5*RHO*L**3 

C 

C       PITCH  HYDRODYNAMIC  COEFHCIENTS 

C 

MQDOT  =-1.680E-02*0.5*RHO*L**5 
MPP  =  5.260E-05*0.5*RHO*L**5 
MPR  =  5.040E-03*0.5*RHO*L**5 
MRR  =-2.860E-03*0.5*RHO*L**5 
MWDOT  =-6.810E-02*0.5*RHO*L**4 
MQ  =-6.860E-02*0.5*RHO*L**4 
MVP     =  1.180E-03*0.5*RHO*L**4 
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MVR     =  1.730E-02*0.5*RHO*L**4 

MW      =  9.860E-02*0.5*RHO*L**3 

MW     =-2.510E-02*0.5*RHO*L**3 

MDS     =-1.113E-02*0.5*RHO*L**3 

MDB     =  1.113E-02*0.5*RHO*L**3 
C 

C       YAW  HYDRODYNAMIC  COEFHCIENTS 
C 

NPDOT  =-3.370E-05*0.5*RHO*L*'^5 

NRDOT  =-3.400E-03*0.5*RHO*L**5 

NPQ     =-2.110E-02*0.5*RHO*L**5 

NQR     =  2.750E-03*0.5*RHO*L**5 

NVDOT  =  1.240E-03*0.5*RHO*L**4 

NP      =-8.405E-04*0.5*RHO*L**4 

NR      =-1.640E-02*0.5*RHO*L**4 

NVQ     =-9.990E-03*0.5*RHO*L**4 

NWP     =-1.750E-02*0.5*RHO*L**4 

NWR     =  7.350E-03*0.5*RHO*L**4 

NV      =-7.420E-03*0.5*RHO*L**3 

NVW     =-2.670E-02*0.5*RHO*L**3 

NDRS    =-1.113E-02*0.5*RHO*L**3 

NDRB    =  +  1.113E-02*0.5*RHO*L**3 
C 

C       VARIABLE  INmUZATION 
C 

RPM=0.0 

BOY = WHGHT+DELB 

XG=XB+XGB 

ZG=ZB+ZGB 

ALPHA=0.0 
C 

C       TDUMMY  IS  DUMMY 

C       TfflS  ROUTINE  EVALUATES  THE  PROBLEM  DEFINING  FUNCTION 
C       INPUT  IS  Y(1),...,Y(N),  AND  PAR=PARAMETER  N0W,F(1),.. 
C       ,F(N)  OF  THE  FUNCTION  OF  RIGHT  HAND  SIDE  OF  F(Y,PAR)  ARE 
C       ENTERED:  (SIX  DEGREE  FREEDOM  EQUATION  OF  MOTION  OF  C 
C       SUBMARINE) 
C 

U=Y(1) 

V=Y(2) 

W=Y(3) 

P=Y(4) 

Q=Y(5) 
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R=Y(6) 
PHI=Y(7) 
THETA=Y(8) 
PSIDOT=Y(9) 
I2=FI2(V,W,Q,R) 
I3=n3(V,W,Q,R) 
I5=FI5(V,W,Q,R) 
I6=FI6(V,W,Q,R) 
SWAY  =12 
HEAVE=I3 
PITCH =15 
YAW =16 
C 

C       PARAMETER  IS  DS  FOR  TfflS  RUN 
C 

DS=PAR 
C 

C       SURGE  FORCE 
C 

F(l)  =  MASS*V*R-MASS*W*Q+MASS*XG*Q**2+MASS*XG*R**2- 
&      MASS*YG*P*Q-MASS*ZG*P*R+XPP*P**2+XQQ*Q**2+XRR*R**2  + 
&      XPR*P*R+XWQ*W*Q+XVP*V*P+XVR*V*R+U*Q* 
&      (XQDS*DS +XQDB*DB)  + 

&      U*R*PCRDRS*DRS+XRDRB*DRB)+XW*V**2+XWW*W**2+U*V* 
&      PCVDRS*DRS+XDRB*DRB)+U*W*(XWDS*DS+XWDB*DB)  + 
&      (XDSDS*DS**2+XDBDB*DB**2+XDRDR*(DRS**2+DRB**2))*U**2- 
&      (WEIGHT-BOY)  *  DSIN(THETA)+XPROP*RPM*RPM-XRES*U*U 
C 

C  SWAY  FORCE 

C 

F(2)  =  -MASS*U*R-MASS*XG*P*Q+MASS*YG*R**2-MASS*ZG*Q*R4- 
&      YPQ*P*Q+YQR*Q*R+YP*U*P+YR*U*R+YVQ*V*Q+YWP*W*P+ 
&      YWR*W*R+YV*U*V+YVW*V*W+YDRS*U**2*DRS+YDRB*U**2* 
&      DRB+  (WEIGHT-BOY)*  DCOS(THETA)*DSIN(PHI)+MASS*W*P+ 
&      MASS*YG*P**2+SWAY 
C 

C  HEAVE  FORCE 

C 

F(3)  =  MASS*U*Q-MASS*V*P-MASS*XG*P*R-MASS*YG*Q*R+ 
&      MASS*ZG*P**2+MASS*ZG*Q**2+ZPP*P**2-hZPR*P*R-h 
&      ZRR*R**2+ZQ*U*Q+ZVP*V*P+ZVR*V*R-I-ZW*U*W-I-ZVV*V**2-H 
&      HEAVE-I-  U**2*(ZDS*DS+ZDB*DB)  +  (WEIGHT-B0Y)* 
&      DCOS(THETA)*  DCOS(PHI) 
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c 

C  ROLL  MOMENT 

C 

F(4)  =  -IZ*Q*R+IY*Q*R-IXY*P*R+IYZ*Q**2-IYZ*R**2+IXZ*P*Q+ 
&      MASS*YG*U*Q-MASS*YG*V*P-MASS*ZG*W*P+KPQ*P*Q+ 
&      KQR*Q*R+  KP*U*P+KR*U*R+KVQ*V*Q+KWP*W*P+KWR* 
&      W*R+KV*U*V+  KVW*V*W+ 

&      (YG*WEIGHT-YB*BOY)*DCOS(THETA)*DCOS(PHI)-(ZG*WEIGHT- 
&      ZB*BOY)*DCOS(THETA)*DSIN(PHI) +MASS*ZG*U*R 
C 

C  PITCH  MOMENT 

C 

F(5)  =  -IX*P*R+IZ*P*R+IXY*Q*R  IYZ*P*Q-IXZ*P**2+IXZ*R**2- 
&      MASS*XG*U*Q+MASS*XG*V*P+MASS*ZG*V*R-MASS*ZG*W*Q+ 
&      MPP*P**2+MPR*P*R+MRR*R**2+MQ*U*Q+MVP*V*P+MVR*V*R 
&      +MW*U*W+MW*V**2+U**2*(MDS*DS+MDB*DB)-(XG*WEIGHT- 
&      XB*BOY)*DCOS(THETA)*DCOS(PHI)- 
&      (ZG*WEIGHT-ZB*BOY)*DSIN(THETA) +PrrCH 
C 

C  YAW  MOMENT 

C 

F(6)  =   IY*P*Q+IX*P*Q+IXY*P**2-IXY*Q**2+IYZ*P*R-IXZ*Q*R- 
&      MASS*XG*U*R+MASS*XG*W*P-MASS*YG*V*R+MASS*YG*W*Q+ 
&      NPQ*P*Q+NQR*Q*R+NP*U*P+NR*U*R+NVQ*V*Q+NWP*W*P+ 
&      NWR*W*R+NV*U*V+  NVW*V*W+NDRS*U**2*DRS+NDRB* 
&      U**2*DRB+(XG*WEIGHT-XB*B0Y)*DC0S(THETA)*DSIN(PHI)  + 
&      (YG*WEIGHT-YB*BOY)*DSIN(THETA) + YAW 
C 

C  KINEMATIC  RELATIONS 

C 

F(7)  =  P+PSIDOT*DSIN(THETA) 
F(8)  =  Q-PSIDOT*DSIN(PHI)*DCOS(THErA) 
F(9)  =  R-PSIDOT*DCOS(PHI)*DCOS(THETA) 
RETURN 
END 

4.       SAMPLE  RUN 

MARS  $  nin  sub2 
WHAT  IS  DRS 
0. 

O.OOOOOOOOOOOOOOOOE+00 

Previous  data  in  DIAGRAM  to  be  purged? 

ENTER  0  FOR  PURGING,  ENTER   1    FOR  APPENDING: 
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0 

THIS  IS  THE  MENU  YOU  CAN  CHOOSE  FROM: 


0  :    ENTERING  STARTING  DATA 

1  :    CONTINUATION:  NO  RESTRICTIONS 

2  :    CONTINUATION:  SETTING  OF  OPTIONS,  THEN  RUN 

3  :    CONTINUATION,  FOLLOWING  PREVIOUS  OPTIONS 
-1    :    ENDING  THIS  SESSION 

ENTER  CHOICE  of  MODUS: 
0 

ENTER  COMPONENT  Y(  1): 
4.2 

ENTER  COMPONENT  Y(  2): 
0 

ENTER  COMPONENT  Y(  3): 
-.7 

ENTER  COMPONENT  Y(  4): 
0 

ENTER  COMPONENT  Y(  5): 
0 

ENTER  COMPONENT  Y(  6): 
0 

ENTER  COMPONENT  Y(  7): 
0 

ENTER  COMPONENT  Y(  8): 
.1745 

ENTER  COMPONENT  Y(  9): 
0 

ENTER  PARAMETER: 
.345 
THIS  IS  THE  MENU  YOU  CAN  CHOOSE  FROM: 


0  :    ENTERING  STARTING  DATA 

1  :    CONTINUATION:  NO  RESTRICTIONS 

2  :    CONTINUATION:  SETTING  OF  OPTIONS,  THEN  RUN 

3  :    CONTINUATION,  FOLLOWING  PREVIOUS  OPTIONS 
-1    :   ENDING  THIS  SESSION 

ENTER  CHOICE  of  MODUS : 
2 

CURRENT  STARTING  DATA  AT  PARAMETER:    0.3450000000000000 
THESE  ARE  THE  CURRENT  OPTIONS: 
No:  1  ,    step  length  :      -0.17406E-01 
No:  2  ,    index  of  fixed  component  :     -1 
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No:  3  ,   level  of  continuation  :     0 

No:  4  ,   maximum  number  of  continuation  steps  :   200 

No:  5  ,   step  bounds:  Predictor: 0.100E+ 00  ,  Param.:0.100E+00  ,  rel. change: 0.10 

No:  6  ,   stability  analysis  (yes=l,  no=0)  :     0 

No:  7  ,   window  :   -0.35000E+00  .It.  parameter  .It.      0.35000E+00 

No:  8  ,   target  point:   no  final  value  is  set. 

ENTER   Number  OF  OPTION  YOU  WANT  TO  CHANGE,   OR 

ENTER  0         TO  RUN  THE  CONTINUATION,    OR 

ENTER   -1       TO  RETURN  TO  MENU: 
1 

ENTER  INITIAL  STEPSIZE  (IN  REAL): 
-.005 

ENTER  Number  OF  OPTION  YOU  WANT  TO  CHANGE,   OR 

ENTER   0        TO  RUN  THE  CONTINUATION,    OR 

ENTER   -1        TO  RETURN  TO  MENU: 
2 

CHOOSE  INDEX  FOR  DSmiAL  STEP: 

ENTER  0  FOR  CONIIN.  WITH  RESPECT  TO  PARAMETER, 
ENTER  K  FOR  FIXING  K-TH  COMPONENT: 
0 

ENTER  Number  OF  OPTION  YOU  WANT  TO  CHANGE,   OR 

ENTER  0         TO  RUN  THE  CONTINUATION,    OR 

ENTER   -1       TO  RETURN  TO  MENU: 
3 

CHOOSE  AMONG  THREE  LEVELS  OF  STEP  CONTROL: 

ENTER  2   FOR  FREELY  VAIOABLE  STEPS, 

1    FOR  VARIABLE  STEP  BUT  FIXED  INDEX, 
0  FOR  BOTH  STEPLENGTH  AND  INDEX  FIXED: 
0 

ENTER  Number  OF  OPTION  YOU  WANT  TO  CHANGE,   OR 

ENTER  0         TO  RUN  THE  CONTINUATION,   OR 

ENTER   -1       TO  RETURN  TO  MENU: 
4 

Enter  maximum  NUMBER  OF  CONTINUATION  STEPS: 
100 

ENTER  Number  OF  OPTION  YOU  WANT  TO  CHANGE,   OR 

ENTER  0         TO  RUN  THE  CONTINUATION,    OR 

ENTER   -1        TO  RETURN  TO  MENU: 
-1 
THIS  IS  THE  MENU  YOU  CAN  CHOOSE  FROM: 


0  :   ENTERING  STARTING  DATA 

1  :    CONTINUATION:  NO  RESTRICTIONS 
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2  :   CONTINUATION:  SETTING  OF  OPTIONS,  THEN  RUN 

3  :    CONTINUATION,  FOLLOWING  PREVIOUS  OPTIONS 
-1    :    ENDING  TfflS  SESSION 

ENTER  CHOICE  of  MODUS: 
2 

CURRENT  STARTING  DATA  AT  PARAMETER:    0.3450000000000000 
THESE  ARE  THE  CURRENT  OPTIONS: 
No:  1  ,    step  length  :      -0.50000E-02 
No:  2  ,   index  of  fixed  component  :      0 
No:  3  ,   level  of  continuation  :     0 
No:  4  ,   maximum  number  of  continuation  steps  :    100 

No:  5  ,    step  bounds:  Predictor: 0.100E+ 00  ,  Param.:0.100E+00  ,  rel.change:0.10 
No:  6  ,    stability  analysis  (yes=l,  no=0)  :     0 
No:  7  ,   window  :   -0.35000E+00  .It.  parameter  .It.      0.35000E+00 
No:  8  ,   target  point:   no  final  value  is  set. 

ENTER  Number  OF  OPTION  YOU  WANT  TO  CHANGE,    OR 

ENTER  0        TO  RUN  THE  CONTINUATION,   OR 

ENTER  -1       TO  RETURN  TO  MENU: 
0 


130 


LIST  OF  REFERENCES 


1.  Clayton,  B.R  and  Bishop,  R.E.D.  Mechanics  of  Marine  Vehicles.  Gulf 
Publishing  Company,  Houston,  1982. 

2.  Roddy,  R.F,  Investigation  of  the  stability  and  control  characteristics  of  several 
configurations  of  the  DARPA  SUBOFF  model  (DRTC  model  5470)  from 
captive-model  experiments,  David  Taylor  Research  Center,  Report  DRTC/SHD- 
1298-08.  1990. 

3.  Falzarano,  J.M.,  Shaw,  S.W,  and  Troesch,  A.W,  "Application  of  Global 
Methods  for  Analyzing  Dynamical  Systems  to  Ship  Rolling  Motion  and 
Capsizing".  Journal  of  Bifurcation  and  Chaos.  2,  1,  1992. 

4.  Taz  Ul  Mulk,  M.  and  Falzarano,  J.,  "Complete  Six  Degrees  of  Freedom 
Nonlinear  Ship  Rolling  Motion".  Ocean  Engineering.  Submitted  for  puUication, 
1993. 

5.  Saipkaya,  T.,  "Brief  Reviews  of  Some  Time-Dependent  Flows",  Journal  of 
Fluids  Engineering.  Vol.  114/283,  September  1992. 

6.  Lugt,  H.J. ,  "Numerical  Modeling  of  Vortex  Flows  in  Ship  Hydrodynamics  -A 
Review",  Proceedings  of  the  Third  International  Conference  on  Numerical  Ship 
Hydrodynamics.  Vol.  1,  1981. 

7.  Booth,  T.B.,  "Stability  of  Buoyant  Underwater  Vehicles  Part  I.  Predominantly 
Forward  Motion".  International  Shipbuilding  Progress.   24,   279,1977. 

8.  Booth,  T.B.,  "Stability  of  Buoyant  Underwater  Vehicles  Part  H.  Near  Vertical 
Ascent",  International  Shipbuilding  Progress.  24,  280,  1977. 

9.  Papoulias,  F.A.,  and  McKinley,  B.D.,  "Inverted  Pendulum  Stabilization  of 
Submarines  in  Free  Positive  Buoyancy  Ascent:,  Journal  of  Ship  Research.  In 
press,  1993. 

10.  Golubitsky,  M.  and  Schaeffer,  D.G.,  Singularities  and  Groups  in  Bifurcation 
Theory  I.  Springer- Verlag,  1985. 


131 


1 1 .  Guckenheimer,  J.  and  Holmes  P. ,  Nonlinear  Oscillations.  Dynamical  Systems. 
and  Bifurcations  of  Vector  Fields.  Applied  Mathematical  Sciences  42,  Springer- 
Verlag,  1983. 

12.  Seydel,  R.,  From  Equilibrium  to  Chaos  Practical  Stability  and  Bifurcation 
Analysis.  Elsevier,  1988. 

13.  Smith,  N.S.,  Crane,  J.W.,  and  Summey,  D.C.,  SDV  Simulator  Hydrodynamic 
Coefficients,  Report  NCSC-TM231-78.  Naval  Coastal  Systems  Center,  1978. 

14.  McKinley  B.D.,  "Dynamic  Stability  of  Positively  Buoyant  Submersibles: 
Vertical  Plane  Solutions",  Master's  Thesis,  Naval  Postgraduate  School, 
Monterey,  California,  December  1991. 

15 .  Seydel,  R. ,  "BIFPACK:  A  Program  Package  for  Calculating  Bifurcations" ,  State 
University  of  New  York  at  Buffalo,  1985. 


132 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 
1 .       Defense  Technical  Information  Center  2 

Cameron  Station 
Alexandria  VA  22304-6145 


Library,  Code  052 
Naval  Postgraduate  School 
Monterey  CA  93943-5002 


3.  Deniz  Harp  Okulu 

Tuzla 
Istanbul,  TURKEY 

4.  Golcuk  Tersanesi  Komutanligi 
Golcuk 

Kocaeli,  TURKEY 

5.  Taskizak  Tersanesi  Komutanligi 
Kasimpasa 

Istanbul,  TURKEY 

6.  Deniz  Kuwetleri  Komutanligi 
Personel  Egitim  Daire  Baskanligi 
Bakanliklar 

Ankara,  TURKEY 

7.  Professor  Fotis  A.  Papoulias,  Code  ME/Pa 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 

Monterey,  CA  93943-5100 

8.  Department  Chairman,  Code  ME/Kk 
Department  Of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey,  CA  93943-5100 


133 


9.  Naval  Engineering  Curricular  Office  (Code  34) 
Naval  Postgraduate  School 

Monterey,  CA  93943-5100 

10.  Ibrahim  AYDIN 

Mobil  Karsisi  Serdar  Apt. 
No:  241/4        Ciftlikkoy 
Yalova         Istanbul/TURKEY 


134 


BUqiEYfCl^oX  LIBRARY 


GAYLORD  S 


